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1  Quantitative  Image  Recovery  From  Measured  Blind 

2  Backscattered  Data  Using  a  Globally 

3  Convergent  Inverse  Method 

4  Andrey  V.  Kuzhuget,  Larisa  Beilina,  Michael  V.  Klibanov,  Anders  Sullivan,  Lam  Nguyen,  and  Michael  A.  Fiddy 


5  Abstract — The  goal  of  this  paper  is  to  introduce  the  application 

6  of  a  globally  convergent  inverse  scattering  algorithm  to  estimate 

7  dielectric  constants  of  targets  using  time-resolved  backscattering 

8  data  collected  by  a  U.S.  Army  Research  Laboratory  for  war  d- 

9  looking  radar.  The  processing  of  the  data  was  conducted  blind,  i.e., 

10  without  any  prior  knowledge  of  the  targets.  The  problem  is  solved 

11  by  formulating  the  scattering  problem  as  a  coefficient  inverse 

12  problem  for  a  hyperbolic  partial  differential  equation.  The  main 

13  new  feature  of  this  algorithm  is  its  rigorously  established  global 

14  convergence  property.  Calculated  values  of  dielectric  constants  are 

15  in  a  good  agreement  with  material  properties,  which  were  revealed 

16  a  posteriori . 

17  Index  Terms — Experimental  data,  inverse  scattering,  quantita- 

18  tive  imaging,  remote  sensing. 

19  I.  Introduction 

20  A  FUNDAMENTAL  problem  in  remote  sensing  is  the 

21  processing  of  scattered  field  data  from  strongly  scattering 

22  penetrable  targets.  Multiple  scattering  renders  this  problem  ex- 

23  tremely  difficult  to  solve,  it  being  ill  conditioned  with  additional 

24  questions  of  uniqueness  and,  the  most  difficult,  nonlinearity  to 

25  contend  with.  In  practice,  limited  noisy  data  typically  require 

26  that  some  physical  models  be  assumed,  from  which  one  hopes 

27  to  extract  meaningful  and  preferably  quantitative  information 

28  about  the  target  in  question.  A  number  of  recent  publications 

29  by  Beilina  and  Klibanov  [3]— [8]  and  by  Klibanov  et  al.  [12], 

30  [14]— [16]  have  led  to  a  new  approach  to  address  this  important 

31  topic.  This  numerical  method  was  originally  developed  for 

32  some  multidimensional  coefficient  inverse  problems  (MCIPs) 

33  for  a  hyperbolic  partial  differential  equation  (PDE)  using  data 
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from  only  a  single  location  of  either  a  point  source  or  from  34 
a  single  direction  of  an  incident  plane  wave.  In  particular,  in  35 
[14],  that  method  was  extended  from  the  3-D  case  to  the  1-D  36 
case.  Thus,  that  1-D  version  of  [14]  is  used  here  to  work  with  37 
the  experimental  data.  The  illuminating  field  is  pulsed  in  time,  38 
and  the  time  history  of  the  backscattering  from  the  illuminated  39 
target  volume  constitutes  the  measured  data  that  are  processed  40 
by  this  algorithm.  The  authors  are  unaware  of  other  groups  41 
working  on  MCIPs  using  data  acquired  from  a  single  source  42 
location.  However,  the  single  measurement  case  is  clearly  the  43 
most  practical  one,  particularly  for  military  applications.  In-  44 
deed,  because  of  many  dangers  on  the  battlefield,  the  number  45 
of  measurements  should  be  minimized.  46 

The  algorithm  in  the  aforementioned  cited  publications  com-  47 
putes  values  for  the  spatial  distribution  of  the  dielectric  con-  48 
stants  of  objects  within  the  target  volume.  It  is  important  to  49 
stress  that  this  algorithm  requires  neither  no  prior  knowledge  50 
of  what  might  exist  in  the  target  volume  nor  a  prior  knowledge  51 
of  a  good  first  guess  about  the  solution.  There  is  a  rigorous  guar-  52 
antee  that  this  algorithm  globally  converges  (see  mathematical  53 
details  in  [7],  [14],  [16],  and  [17]).  Because  of  the  global  con-  54 
vergence  property,  estimates  of  spatially  distributed  dielectric  55 
constants  are  reliable  and  systematically  improve  with  more  56 
measured  and  less  noisy  data.  The  theory  of  the  aforementioned  57 
cited  publications  rigorously  guarantees  that  this  numerical  58 
method  delivers  a  good  approximation  to  the  exact  solution  59 
of  an  MCIP  without  any  a  priori  information  about  a  small  60 
neighborhood  of  the  exact  solution  as  long  as  iterations  start  61 
from  the  so-called  “first  tail  function”  V${x),  which  can  be  62 
easily  computed  using  available  boundary  measurements  (see  63 
(2.27)-(2.29)  in  Section  II-C).  In  addition,  it  is  in  this  sense  64 
that  we  use  the  term  “global  convergence"  of  the  algorithm.  65 
The  common  perception  of  the  term  "global  convergence"  is  66 
that  one  can  start  from  any  point  and  still  get  the  solution,  but  67 
we  stress  that  we  actually  start  not  from  any  point  but  rather  68 
from  the  function  Vo(x ),  which  can  be  easily  computed  from  69 
the  boundary  data  (see  (2.27)-(2.29)  in  Section  II-C).  70 

It  is  well  known  that  least  squares  functionals  for  MCIPs  71 
suffer  from  multiple  local  minima  and  ravines.  Hence,  local  72 
convergence  of  numerical  methods  to  incorrect  estimates  will  73 
occur  unless  an  initial  guess  that  is  close  to  the  true  solution  is  74 
used.  Such  a  guess  is  rarely  available  in  most  applications.  In  75 
contrast,  our  algorithm  does  not  use  a  least  squares  functional,  76 
and  hence,  it  is  free  from  the  problem  of  local  minima.  Instead,  77 
this  algorithm  relies  on  the  structure  of  the  differential  operator  78 
of  the  wave-like  PDE.  79 
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80  Prior  to  the  work  reported  here,  a  major  focus  by  the 

81  U.S.  Army  Research  Laboratory  (ARL)  had  been  on  the  de- 

82  velopment  of  image  processing  techniques  [19]  that  would 

83  improve  radar  images,  which  is  through  postprocessing  tech- 

84  niques  rather  than  through  the  application  of  inverse  scattering 

85  methods.  By  incorporating  more  physics  of  the  target- wave 

86  electromagnetic  response  into  the  data  processing,  one  can 

87  greatly  improve  target  detection  and  identification.  Present  data 

88  processing  provides  an  electromagnetic  field  brightness  or  an 

89  intensity  map  of  the  target  volume,  which  need  not  relate 

90  in  a  simple  fashion  to  the  scattering  structures  themselves. 

91  Our  method  estimates  dielectric  constants  of  targets,  which 

92  obviously  adds  an  important  new  dimension  to  the  interpre- 

93  tation  of  data  acquired  by  the  radar  system  since  this  allows 

94  specific  bounds  on  the  dielectric  properties  of  a  feature  in 

95  the  target  volume,  which  can  help  identify  its  likely  material 

96  properties.  Since  no  prior  knowledge  is  required,  the  measured 

97  data  were  processed  by  Kuzhuget,  Beilina,  Klibanov,  and  Fiddy 

98  in  the  most  challenging  scenario,  i.e.,  without  any  knowledge 

99  of  the  actual  target  structures  and  their  dielectric  properties. 

100  Once  this  had  been  done,  Sullivan  and  Nguyen  compared  a 

101  posteriori  the  image  estimates  with  the  actually  known  material 

102  characteristics. 

103  We  draw  attention  to  the  fact  that  this  algorithm  has  been 

104  used  with  forward-scattered  data  from  experiments.  These 

105  results  were  previously  reported,  which  are  also  in  a  blind 

106  experiment  (see  [12,  Tables  5  and  6]  and  [7,  Tables  5.5  and 

107  5.6]).  In  this  case,  the  images  in  [12]  were  further  improved 

108  and  presented  in  a  follow-up  publication  [6]  using  the  adaptivity 

109  technique  of  [1],  [2],  [4],  [5],  and  [7]. 

110  In  Section  II,  we  outline  the  basic  steps  in  the  underlying 

111  theory  upon  which  the  new  algorithm  is  based.  In  Section  III, 

112  we  formulate  the  global  convergence  theorem.  In  Section  IV, 

113  we  outline  results  obtained  using  time-resolved  backscatter 

114  electric  field  measurements  collected  in  the  field.  Measure- 

115  ments  were  carried  out  by  a  forward-looking  radar  system  built 

116  and  operated  by  the  ARL.  The  data  were  noisy  and  limited,  and 

117  the  target  volumes  included  miscellaneous  sources  of  clutter. 

118  The  purpose  of  this  particular  radar  system  is  to  detect  and 

119  possibly  identify  shallow  explosive-like  targets. 


120 


II.  Theoretical  Background 


We  assume  that  the  constitutive  parameter  of  interest,  i.e.,  135 
mapping  the  target  volume,  is  a  relative  permittivity  er(x).  In  136 
other  words,  we  ignore  magnetic  effects  in  this  paper.  We  also  137 
assume  for  convenience  that  er(x)  =  1  outside  of  the  target  138 
volume,  which  is  x  £  (0, 1)  in  our  case.  We  assume  that  the  139 
source  xo  <  0  lies  outside  of  the  target  volume.  We  can  write  140 
the  forward  scattering  problem  as  141 

Er{pc)'Ujtt  — xx •>  X  £  M  (2.1) 

?i(x,0)=0,  ut(x,  0)  =  S(x  —  xo).  (2.2) 

The  subscripts  in  (2.1)  indicate  the  number  of  partial  derivatives  142 
with  respect  to  the  variable  indicated.  The  coefficient  inverse  143 
problem  (CIP)  is  to  recover  er(x ),  assuming  that  the  initial  144 
illuminating  pulse  is  known  and  that  we  measure  the  function  145 
g(t ),  i.e.,  146 


u(0,t)  =  g(t) 


(2.3) 


121  A.  Integral  Differential  Equation 

122  Since  we  were  given  only  one  time-resolved  experimental 

123  curve  per  target,  we  had  no  choice  but  to  use  a  1-D  mathemati- 

124  cal  model,  although  the  reality  is  3-D  (see  Section  III  for  some 

125  details  about  the  data  collection).  In  addition,  since  only  one 

126  component  of  the  electric  wave  field  was  both  transmitted  and 

127  measured,  we  model  the  scattering  process  with  one  wave-like 

128  PDE  rather  than  using  complete  Maxwell  equations.  We  stress 

129  that  the  method  is  designed  for  use  with  3-D  problems,  and 

130  one  would  normally  collect  data  with  co  polarization  and  cross 

131  polarization  in  order  to  capture  all  of  the  pertinent  information 

132  about  the  target.  Here,  we  simply  wish  to  show  the  steps 

133  employed  by  the  method  and  demonstrate  their  quantitative 

134  reconstruction  accuracy  given  noisy  measured  data. 


for  sufficiently  large  times  t  that  all  multiple  scattering  events  147 
within  the  target  volume,  which  can  produce  a  measurable  148 
signal  at  the  detector,  do  so.  Practically,  we  gate  the  radiation  149 
source  in  time;  and  since  the  Laplace  transform  (LT),  i.e.,  150 
w(x,s),  is  used  to  solve  this  CIP,  the  decay  e~st,s  >  0  of  151 
the  LT  kernel  further  limits  the  duration  of  the  measured  time  152 
history.  It  is  worth  pointing  out  that,  more  typically,  scattering  153 
data  would  be  measured  at  different  scattering  angles  for  fixed  154 
frequency  illumination  at  various  incident  angles.  One  can  155 
easily  appreciate  that  this  leads  to  the  acquisition  of  Fourier  156 
information  about  the  target  or  the  secondary  source  function,  157 
depending  upon  the  extent  of  the  multiple  scattering;  and  once  158 
one  has  sufficient  data,  a  reasonable  estimate  of  the  target  159 
properties  becomes  possible.  By  taking  measurements  in  the  160 
time  domain,  one  can  see  that  this  is  essentially  simultane-  161 
ously  gathering  information  in  a  transform  space  from  many  162 
illumination  frequencies.  The  Laplace  and  Fourier  transforms  163 
provide  complimentary  representations  of  the  target  in  terms  of  164 
moments  or  modes,  respectively.  165 

The  LT  is  166 

03 

M*,s)  =  ju(x,t}e-<,t:=Cu,  »>s=conS,.>0  (2.4) 

0 

and  we  assume  that  the  so-called  pseudofrequency  s  >  167 
s(er(x))  :=  s  is  sufficiently  large.  This  gives  [7]  168 

wxx  ~  s2er(x)w  =  —  S(x  —  xq),  x  £  M  (2.5) 


lim  w(x,  s )  =  0. 

|^-0O3 


(2.6) 


Let 


169 


w(  0,  s)  =  ip(s)  =  Cg 


(2.7) 


be  the  LT  of  the  measured  function  g(t)  in  (2.3).  Since  er(x)  =  170 
1  for  x  <  0,  then,  using  (2.5)  and  (2.6),  one  can  prove  that,  in  171 
addition  to  the  function  rc(0,  s)  in  (2.7),  the  function  wx( 0,  s)  172 
is  also  known  as  (see  [17])  173 


Wx( 0,  s)  =  Sip(s)  -  exp(sxo). 


(2.8) 
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174  Let  wo  (x,  s)  be  the  solution  of  the  problem  in  (2.5)  and  (2.6) 

175  for  the  case  of  the  uniform  background  er(x)  =  1.  Then 

exp  (— s\x  —  #o|) 


WQ(X,S)  = 


2s 


r(x,s)  =  4ln  (  —  (x,s)  )  . 

Sz  \Wo 


qxx  +  2s2qxrx  +  2srl~2sqx-2rx=0, 


V (x,  s )  :=  V (x)  =  r(x,  s)  =  4  ln  ( 

s  \Wo(x,S) 


Qxx  2s  qx 

—  2 sqx  +  2 


ir  +  2s 


J  qx(x,r)d 

S 

s 

J  qx(x,T)dT 

S 

S 

+  2s2qxVx  -  4sVx  J  qx(x,  r)d 

+  2s{Vxf-2Vx  =0, 
x  G  (0, 1);  s  G  [s,  s] 
q(0,s)  =  ipo(s),  qx(0,s)  =  i(>1(s) 
qx( l,s)sO,  s  G  [s,  s] 


J  qX{x,T)d7 


(2.9) 


176  When  implementing  the  algorithm,  given  the  assumption  of  a 

177  uniform  normalized  er(x)  =  1  outside  of  the  target  volume,  we 

178  consider  the  function 


(2.10) 


179  Since  the  source  xo  <  0,  then  the  function  r(x,  s )  is  the  solution 

180  of  the  following  equation  in  the  interval  (0,  1): 

rXx  +  s2r2x  -  2srx  =  er(x)  -  1,  iG(0,l).  (2.11) 

181  In  addition,  by  (2.7)  and  (2.8) 

r(0,s)  =ip0(s),  rx(0,  s)  =  ipi(s)  (2.12) 

ln</?(s)  -  ln(2s)  x0 

lPo{s)  ~ - ~2 - ^  GT 

sz  s 

2  esx° 

<Pi(*)= - 2-TT-  (2'13) 

S  s2ip(s) 

182  The  idea  now  is  to  eliminate  the  unknown  coefficient  er(x) 

183  from  (2.11)  via  differentiation  with  respect  to  pseudofre- 

184  quency  s.  Differentiating  (2.11)  with  respect  to  s  and  denoting 

185  q(x,  s )  =  dsr(x ,  s),  we  obtain 


x  G  (0, 1).  (2.14) 


186  We  now  need  to  express  in  (2. 14)  the  function  r  via  the  function 

187  q.  We  have 

s 

>■(*.»)  =  -  j  q(x,T)dT  +  V(x,H)  (2.15, 


188  where  V{x)  :=  V  (x,  s)  is  referred  to  as  the  tail  function,  which 

189  is  small  in  practice  for  large  positive  s.  Here,  the  truncation 

190  pseudofrequency  s  serves  as  a  regularization  parameter.  The 

191  exact  formula  for  V  (x)  is 


Fig.  1.  Flowchart  of  the  algorithm. 

where  functions  ^o( <5 )  =  <Po(s)  and^i(s)  =  ipi(s)  are  derived  194 
from  (2.13).  The  condition  qx(l,  s)  =  0  can  be  easily  derived  195 
from  (2.6)  since  er{pc)  =  1  outside  of  the  interval  (0,  1).  196 

In  (2.17)  and  (2.18),  both  functions  q(x,s)  and  V(x)  are  197 
unknown.  The  reason  why  we  can  approximate  both  of  them  198 
is  that  we  find  updates  for  q(x,s)  via  inner  iterations  exploring  199 
(2.17)  and  (2.18)  inside  of  the  interval  (0,  1).  At  the  same  time,  200 
we  update  the  tail  function  V  (x)  via  outer  iterations  exploring  201 
the  entire  real  line  M.  In  short,  given  an  approximation  for  V (x),  202 
the  algorithm  updates  q  and  then  updated  for  er(x).  Next,  the  203 
forward  problem  in  (2.5)  and  (2.6)  is  solved  for  the  function  204 
w(x,s)  for  s  =  s.  Next,  the  tail  function  V (x)  is  updated  using  205 
(2.16).  This  might  seem  reminiscent  of  the  steps  in  algorithms  206 
such  as  the  modified  gradient  inverse  scattering  technique  [20] ;  207 
but  we  emphasize  that,  unlike  our  case,  such  methods  have  no  208 
global  convergence  properties.  209 


B.  Iterative  Process 


210 


(2.16) 


192  Substituting  (2.15)  in  (2.14),  we  obtain  the  following  nonlinear 

193  integral  differential  equation: 

-i  2 


We  now  outline  the  formulation  of  our  algorithm  and  the  211 
iterative  process  (see  details  in  [7],  [14],  [16],  and  [17];  see  212 
Fig.  1).  Unlike  computationally  simulated  data  in  [14],  we  213  AQ2 
do  not  use  prior  knowledge  of  the  function  g(l,s)  on  the  214 
transmitted  edge  since  this  function  is  unknown  to  us.  We  have  215 
observed  in  our  computational  experiments  that  the  knowledge  216 
of  q(l,s)  only  affects  the  accuracy  of  the  calculation  of  the  217 
location  of  the  target,  but  it  does  not  affect  the  accuracy  of  the  218 
computed  target/background  contrast.  Here,  we  are  interested  219 
only  in  that  contrast  (see  Section  III).  Since  er(x)  =  1  for  x  >  1  220 
and  xo  <  0,  then  one  can  easily  derive  from  equations  (2.5),  221 
(2.9),  and  (2.10)  that  dxq(  1,  s )  =  0.  222 

Consider  a  partition  of  the  interval  [s,  s]  into  N  small  subin-  223 
tervals  with  the  small  grid  step  size  h  >  0  and  assume  that  the  224 
function  q(x,  s)  is  piecewise  constant  with  respect  to  s.  Thus  225 


S  =  S tv  <  Sjv-1  <  •  •  •  <  So  =  5, 
q(x,s)  =qn(x),  for  s  G  (sn,  sn_i]. 


Si- 1  —  Si  =  h 

(2.19) 


(2.17) 

(2.18) 


For  each  subinterval  (sn,  sn-i]  we  obtain  a  differential  equation  226 
for  the  function  qn(x).  We  assign,  for  convenience  of  notations,  227 
qo  :=  0.  Following  the  aforementioned  idea  of  a  combination  of  228 
inner  and  outer  iterations,  we  perform  for  each  n  inner  iterations  229 


230  with  respect  to  the  tail  function.  This  way,  we  obtain  functions 

231  </„,/,  and  Vn,k.  The  equation  for  the  pair  14, fc)  is 


n— 1 


Qn,k  1  A\ ,n/z  ^  ]  Qj  AijnVn,k  2^2, n  I  Qn,k 


3= 0 
^  n—1 


n— 1 


n—1 


~A2,nh2  {'El]  +2h^j+2A2tnV^k\h^j 


vi=° 
\  2 


i=o 


i=o 


V'  , 

i  v  n^ki 

xe(0,i) 

O 

Jl^ 

t—i 

3 

II 

1 

sn—  1 
/» 

(2.21) 


232  Here,  A\,n  and  A2,n  are  certain  numbers,  whose  exact  expres- 

233  sions  are  given  in  [3]  and  [7]. 

234  The  choice  of  the  first  tail  function  Vo(x)  is  described  in 

235  Section  II-C.  Let  n  >  1.  Suppose  that,  for  j  =  0, . . .  n  —  1, 

236  functions  qj{x)  and  Vj(x)  are  already  constructed.  We  now 

237  need  to  construct  functions  qn^  and  Vnjk  for  k  =  1, . . . ,  m. 

238  We  set  Vn,i(x)  :=  Vn-\{x).  Next,  using  the  quasi-reversibility 

239  method  (QRM)  (see  Section  II-C),  we  approximately  solve 

240  (2.20)  for  k  —  1  with  overdetermined  boundary  conditions  in 

241  (2.21)  and  find  the  function  qn^.  Next,  we  find  the  approxima- 

242  tion  for  the  unknown  coefficient  er{x)  via  the  following 

243  two  formulas: 

n—l 

rn, lW  =  -  hqn, i  -  h^qj  +  14, l,  x  G  [0,1]  (2.22) 

J=0 

4n,1)(X)  =  i  +  r"^)  +  4  [4,i  60] 2 

-2snr'nA(x),  x  e  [0, 1].  (2.23) 

244  Next,  we  solve  the  forward  problem  in  (2.5)  and  (2.6)  with 

245  er(x)  :=  £rn,1)(^),  s  :=s  and  find  the  function  wn^(x,s) 

246  this  way.  After  this,  we  update  the  tail  via  the  formula  in  (2.16), 

247  in  which  w(x,  s)  :=  wn^(x,  s).  This  way,  we  obtain  a  new  tail 

248  Vn^{x).  Similarly,  we  continue  iterating  with  respect  to  tails  m 

249  times.  Next,  we  set 

qn(x)  ■■=  qn,m(x),  Vn(x)  :=  Vn,m(x),  e^\x :)  := 

250  replace  n  with  n  +  1  and  repeat  this  process.  We  continue  this 

251  process  until  [15] 


either 


£(n)  _  ^O-1) 


MOT) 


<  10 


-5 


or  ||VJa(gnjfe)||L2(0jl)  >  105  (2.24) 


252  where  the  functional  Ja{qn,k )  is  defined  in  Section  II-C.  Here, 

253  the  norm  in  the  space  L2(0, 1)  is  understood  in  the  discrete 

254  sense.  In  the  case  when  the  second  inequality  in  (2.24)  is 

255  satisfied,  we  stop  at  the  previous  iteration,  taking  e (x)  as 

256  our  solution.  If  neither  of  two  conditions  in  (2.24)  is  not  reached 

257  at  n  :=  N,  then  we  repeat  the  aforementioned  sweep  over  the 

258  interval  [s ,  s] ,  taking  the  pair  (q]y(x),  Vn(x))  as  the  new  pair 

259  (qo(x),  V& (a)).  Usually,  at  least  one  of  the  conditions  in  (2.24) 

260  is  reached  either  on  the  third  or  on  the  fourth  sweep,  and  the 
AQ3  261  process  stops  then. 
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C.  Computing  Functions  qn,k{x)  and  Vq(x)  262 

At  first  glance,  it  seems  that,  for  a  given  tail  function  Vn^  (x),  263 
the  function  qn,k(x)  can  be  computed  as  the  solution  of  a  264 
conventional  boundary  value  problem  for  the  ordinary  differ-  265 
ential  equation  in  (2.20)  with  any  two  out  of  three  boundary  266 
conditions  in  (2.21).  However,  attempts  to  do  so  led  to  poor  267 
quality  images  (see  [14,  Remark  3.1]).  At  the  same  time,  the  268 
QRM  has  resulted  in  accurate  solutions  both  in  [14]  and  in  Test  269 
1  for  synthetic  data  (see  succeeding  discussion).  The  QRM  is  270 
well  designed  to  compute  least  squares  solutions  of  PDEs  with  271 
overdetermined  boundary  conditions,  such  as,  e.g.,  the  problem  272 
in  (2.20)  and  (2.21).  We  refer  to  [18]  for  the  originating  work  273 
about  the  QRM  and  to  [7],  [9],  [13],  [15],  and  [16]  for  some  274 
follow-up  publications.  275 

Let  L(qn,k)(x)  and  Pn^{x)  be  left-  and  right-hand  sides  of  276 
(2.20),  respectively.  In  our  numerical  studies,  L(qn^)(x)  and  277 
PU:k{x)  are  written  in  the  form  of  finite  differences.  Let  a  G  278 
(0, 1)  be  the  regularization  parameter.  The  QRM  minimizes  the  279 
following  Tikhonov  regularization  functional:  280 

Jaijln^k)  =  ||  Ln,k  ( qn,k )  Pn,k  IIl2(0,1)  ^ll^n,fc  II  JT2(0,1)  (2-25) 

subject  to  boundary  conditions  in  (2.21).  Here,  again  norms  281 
in  L2(0, 1)  and  in  the  Sobolev  space  H2( 0, 1)  are  understood  282 
in  the  discrete  sense.  The  functional  Ja(qn,k )  in  (2.25)  is  283 
quadratic.  Using  this  fact  and  the  tool  of  Carleman  estimates,  it  284 
can  be  shown  that  Ja(qn,k )  has  a  unique  global  minimum  and  285 
no  local  minima  [14],  [15],  [17].  We  find  that  global  minimum  286 
via  the  conjugate  gradient  method,  minimizing  with  respect  to  287 
the  values  of  the  function  qn^  at  grid  points.  We  have  used  288 
100  grid  points  in  the  interval  (0,  1).  The  step  size  in  the  s-  289 
direction  was  h  =  0.5.  The  8-interval  was  [s,  s]  =  [3, 12].  For  290 
each  n  =  1, . . . ,  N,  we  take  functions  qn^  for  k  =  1, ...  m,  291 
and  we  typically  choose  m  =  10.  The  reason  for  the  choice  292 
of  m—  10  is  that  numerical  experience  has  shown  that,  for  293 
each  of  the  n ,  tails  stabilize  at  k  «  10.  As  to  the  regularization  294 
parameter  a ,  we  have  found,  when  testing  synthetic  data,  that  295 
ot  =  0.04  is  the  optimal  one,  and  we  take  it  in  our  computations.  296 

We  note  that  we  determined  the  regularization  parameter  297 
when  testing  simulated  data.  These  data  were  for  the  target  298 
depicted  in  Fig.  7(a),  for  which  we  varied  the  regularization  299 
parameter  between  0.03  and  0.05.  The  resulting  images  for  300 
these  data  showed  only  an  insignificant  change.  We  also  var-  301 
ied  the  regularization  parameter  between  0.03  and  0.05  for  302 
the  experimental  data.  Again,  we  only  observed  insignificant  303 
changes,  which  lead  us  to  select  the  average  value  of  0.04.  304 
Although  the  regularization  theory  states  that  the  regularization  305 
parameter  should  depend  on  the  noise  level  in  the  data  [23],  we  306 
do  not  actually  know  the  noise  level  for  our  data.  Further,  for  307 
nonlinear  problems  (as  we  have),  this  dependence  is  claimed  308 
by  regularization  theory  only  for  the  limiting  case  of  a  relatively  309 
small  level  of  noise,  which  is  not  our  case.  In  our  computations  310 
using  measured  data,  one  works  with  some  level  of  noise,  which  311 
is  not  likely  to  be  small  and  is  unknown.  Therefore,  in  practice,  312 
when  applying  this  algorithm  to  experimental  data,  we  were  313 
guided  by  results  from  simulations  to  choose  a  value  for  the  314 
regularization  parameter.  If  we  had  prior  knowledge  about  some  315 
objects  in  the  target  volume,  then  we  would  choose  the  optimal  316 
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317  regularization  parameter  for  that  object.  Because  we  processed 

318  the  data  without  any  prior  knowledge  whatsoever  about  the 

319  objects,  we  chose  the  regularization  parameter  based  on  the 

320  simulated  data  processing,  and  fortunately,  our  answers  for  five 

321  out  of  five  targets  were  well  within  tabulated  limits. 

322  We  now  describe  an  important  step  in  choosing  the  first 

323  tail  function  Vq(x).  To  choose  it,  we  consider  the  asymptotic 

324  behavior  of  the  function  V(x,s)  in  (2.16)  with  respect  to  the 

325  truncation  pseudofrequency  s  -A  oo.  This  behavior  is  [14],  [17] 


V(x,s)  =  V-^ 

s 


0IP 


oo. 


326  We  truncate  the  term  0(1/ s2),  which  is  somewhat  similar  with 

327  the  defining  of  geometrical  optics  as  a  high-frequency  approx- 

328  imation  of  the  solution  of  the  Helmholtz  equation.  Hence,  we 

329  take 


V(x,s) 


Po(x) 


330  Since  q  =  dsr  and  V (x,  s )  =  r(x,  s ),  then 

Po(x) 


q(x,s)  = 


s2 


(2.26) 


331  Hence,  setting  in  (2.17)  s  :=  s  and  using  (2.26),  we  obtain  the 

332  following  approximate  equation  for  the  function  po(x): 


d 2 

^»(*)  =  o, 


X  £  (0, 1). 


(2.27) 


V  o(-r)  := 


(2.29) 


IV’o(s)  -  $S(S) I  <  <7,  IV’i(s)  -  V’iOOI  < 


for  s  E  Is,  s } 


346  where  functions  'ipo(s)  and  Vh  (s)  depend  on  the  function  g(t)  in 

347  (2.3)  via  (2.7),  (2.13)  and  (2.18);  and  functions  'iPq(s)  and  'ipl(s) 

348  depend  on  the  noiseless  data  g*  ( t )  in  the  same  way.  Let  h  E 

349  (0, 1)  be  the  grid  step  size  in  the  5-direction  in  (2.19);  let  yj a  = 

350  cr  and  h  =  max(cr,  h).  Let  Q  be  the  total  number  of  functions 

(n  k) 

351  Sr  ’  computed  in  the  aforementioned  algorithm.  Then,  there 

352  exists  a  constant  D  =  D(xo,  d,  s)  >  1  such  that,  if  the  numbers 

353  a  and  h  are  so  small,  that 

^  <  JJ2Q+2  (2.30) 


Fig.  2.  (a)  Schematic  diagram  of  the  forward-looking  radar  system  illumi¬ 

nating  a  dielectric  target,  (b)  Typical  measured  time  history  of  the  backscatter 
held,  (c)  Composite  of  unprocessed  returns  highlighting  the  dielectric  target 
(indicated  by  the  red  triangle),  (d)  Downrange  cut  of  the  permittivity  profile, 
which  the  new  algorithm  will  generate. 


then  the  following  estimate  is  valid: 

r(n,/c)  _ 


354 


L2(  0,1) 


<  hu 


(2.31) 


_(n,k) 


333  Boundary  conditions  for  po(x)  can  be  easily  derived  from 

334  (2.18)  and  (2.26)  as 

Po(0)  =  ~s2Ms),  Po(0)  =  ~s2Ms),  PbO)  =  0-  (2.28) 

335  We  find  an  approximate  solution  po,appr(x)  of  the  problem  in 

336  (2.27)  and  (2.28)  via  the  QRM,  similarly  with  the  aforemen- 

337  tioned  equation.  Next,  we  set  for  the  first  tail  function,  i.e., 

P0,appr  (x) 


338  A  simplified  formal  statement  of  the  global  convergence 

339  theorem  is  as  follows  (see  [7,  Th.  6.1]  for  more  details  and 

340  [7,  Th.  6.7]  for  the  3-D  case). 

341  Theorem  1:  Let  the  function  e*r(x)  be  the  exact  solution  of 

342  our  CIP  for  the  noiseless  data  g*(t)  in  (2.3).  Fix  the  truncation 

343  pseudofrequency  s  >  1.  Let  the  first  tail  function  Vq(x)  be 

344  defined  via  (2.27)-(2.29).  Let  cr  E  (0, 1)  be  the  level  of  the  error 

345  in  the  boundary  data,  i.e., 


where  the  number  uj  E  (0, 1)  is  independent  of  n,  k ,  h,  £r  ’"A  355 
and  e*.  356 

Therefore,  Theorem  1  guarantees  that,  if  the  total  number  357 

(n  k) 

Q  of  computed  functions  ;  is  fixed  and  error  parameters  358 
cr,  h  are  sufficiently  small,  then  obtained  iterative  solutions  359 
£rn’^  ( x )  are  sufficiently  close  to  the  exact  solution  e* ,;  and  this  360 
closeness  is  defined  by  the  error  parameters.  Therefore,  the  total  361 
number  of  iterations  Q  can  be  considered  as  the  regularization  362 
parameter  of  our  process,  which  is  the  additional  regularization  363 
parameter  to  the  number  s.  The  combination  of  inequalities  364 
in  (2.30)  and  (2.31)  has  a  direct  analog  in  the  inequality  in  365 
[11,  Lemma  6.2,  p.  156]  for  classical  Landweber  iterations,  366 
which  are  defined  for  a  substantially  different  ill-posed  prob-  367 
lem.  As  to  the  total  number  of  iterations  Q  being  a  regulariza-  368 
tion  parameter  here,  there  is  no  surprise  in  this.  Indeed,  it  is  369 
stated  on  [1 1,  p.  157]  that  the  number  of  iterations  can  serve  as  370 
a  regularization  parameter  for  an  ill-posed  problem.  371 


III.  Imaging  Results 


372 


The  schematic  of  the  data  collection  by  the  forward-looking  373 
radar  is  shown  in  Fig.  2(a).  Time-resolved  electromagnetic  374 
pulses  are  emitted  by  two  sources  installed  on  the  radar.  Only  375 
one  component  of  the  electric  field  is  both  transmitted  and  376 
measured  in  the  backscatter  direction.  The  data  are  collected  377 
by  sixteen  detectors  with  the  step  size  in  time  of  0.133  ns.  378 
Data  from  shallow  targets  placed  both  below  and  above  the  379 
ground  were  provided.  The  only  piece  of  information  provided  380 
by  the  ARL  team  (Sullivan  and  Nguyen)  to  Kuzhuget,  Beilina,  381 
Klibanov,  and  Fiddy  was  whether  the  target  was  located  above  382 
the  ground  or  was  buried.  The  depth  of  the  upper  surface  of  a  383 
buried  target  was  a  few  centimeters.  GPS  was  used  to  provide  384 
the  distance  between  the  radar  and  a  point  on  the  ground,  which  385 
is  located  above  that  target  to  within  a  few  centimeters  error.  386 
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Fig.  3.  (a)  Scattered  field  from  a  metallic  target,  (b)  Scattered  field  from  a  high-permittivity  target  with  the  same  shape  (er  (target)  =  10).  Note  the  similarity 

between  the  backscatter  electric  fields  in  cases  (a)  and  (b). 


387  The  time-resolved  voltages  induced  by  the  backreflected  signals 

388  were  integrated  over  the  radar  to  target  distances  ranging  from  8 

389  to  20  m,  and  they  were  also  averaged  with  respect  to  both  source 

390  positions  and  with  respect  to  the  output  of  the  16  detectors. 

391  Since  we  can  assume  here  that  the  radar/target  distance  was 

392  known,  then  it  was  also  approximately  known  which  part  of  the 

393  measured  time-resolved  signal  would  correspond  to  scattering 

394  events  from  that  target  (see  Fig.  2).  Despite  the  presence  of 

395  clutter,  a  single  time-dependent  curve  is  extracted  from  the 

396  measured  return  time  histories,  as  illustrated  in  Fig.  2(b).  This 

397  is  the  form  of  the  data  that  have  been  processed  in  each  of 

398  the  five  measured  data  sets  generated  by  the  ARL.  A  typical 

399  plot  of  returns  without  applying  the  inverse  algorithm  is  shown 

400  in  Fig.  2(c),  where  the  triangle  denotes  a  possible  target  of 

401  interest  among  the  clutter  from  the  backscatter  generated  from 

402  the  volume  of  the  region  illuminated  by  the  radar  in  Fig.  2(a). 

403  We  process  a  set  of  averaged  time  histories  like  those  shown  in 

404  Fig.  2(b)  to  create  a  down-range  cut  of  the  permittivity  profile, 

405  as  indicated  in  Fig.  2(d). 

406  Our  objective  was  to  calculate  ratios 


R  = 


er  (target) 
n  (background) 


(3.1) 


illustrates  this.  We  use  the  upper  bound  sr  (target)  =  30  for  425 
the  metallic  targets  because  our  calculations  show  that  LT  in  426 
(2.7),  from  the  response  function  g(t),  almost  coincides  for  427 
er  (target)  >  30.  428 

In  both  cases  of  a  metal  structure  and  a  high-permittivity  429 
structure,  one  can  expect  enhanced  backscatter  if  the  incident  430 
pulse  includes  frequencies  that  correspond  to  a  normal  mode  of  431 
the  target.  Hence,  we  assign  432 


10  <  er  (metallic  target)  <  30. 


(3.2) 


407  of  dielectric  constants.  If  the  er  (background)  is  known,  then  it 

408  is  trivial  to  deduce  er  (target) .  Clearly,  for  a  target  located  above 

409  the  ground,  er  (background)  =  1.  In  general,  we  would  expect 

410  the  target  volume  to  contain  many  inhomogeneities  with  spa- 

411  tially  varying  sr(x).  A  weighted  average  of  dielectric  constants 

412  of  these  constituent  materials  will  be  found  over  the  volume 

413  spatial  resolution  cell  that  corresponds  to  the  particular  data 

414  acquisition  configuration.  In  the  examples  presented  here,  we 

415  show  results  obtained  from  just  one  time-history  curve  for  each 

416  target,  corresponding  to  one  polarization  component  of  the  in- 

417  cident  electromagnetic  field  and  backscatter  data  measured  and 

418  averaged  over  all  16  receiver  locations.  Clearly,  this  severely 

419  limits  the  transverse  resolution  but  improves  the  signal-to-noise 

420  ratio  for  1-D  imaging  in  the  propagation  direction.  The  model 

421  is  further  simplified  by  using  the  1-D  CIP  employing  only 

422  one  hyperbolic  PDE.  Consequently,  the  interpretation  of  the 

423  backscattering  radiation  will  assign  a  high-permittivity  value 

424  to  metal  structures.  A  comparison  between  Fig.  3(a)  and  (b) 


We  call  (3.2)  the  appearing  dielectric  constant  of  metallic  tar¬ 
gets.  In  other  words,  we  consider  in  (3.2)  that  regions  appearing 
to  have  a  high  dielectric  constant  could  also  be  metallic  targets. 

To  appreciate  the  kind  of  backscatter  data  and  image  recov¬ 
ery  expected  from  a  simple  dielectric  block,  a  1-D  example 
illustrated  in  Fig.  3  was  investigated.  Computations  in  this 
example  were  performed  using  the  software  package  WavES 
[24].  The  permittivity  profile,  i.e.,  sr  (target)  =  4,  is  shown  in 
Fig.  4(a);  and  the  computed  function  u(0,t)  =  g(t)  for  0  < 
t  <  3  is  shown  in  Fig.  4(b)  [see  (2.3)  for  g(t)].  We  assume 
temporal  units  here  for  which  at  t  =  3,  a  distance  of  x  =  3 
units  is  traversed;  the  source  is  at  xq  =  —  1,  and  the  block’s 
front  face  is  at  x  =  0.2.  Since  the  block  is  0.2  units  wide,  g(t) 
represents  the  backscatter  return  from  the  front  and  back  face  of 
the  block.  The  reason  why,  in  Fig.  4(b),  g(t )  =  0  for  t  <  1  and 
g(t)  =  1/2  for  1  <t<  1.4  is  that  the  solution  of  the  problem  in 
(2.1)  and  (2.2)  for  er(x)  =  1  is  uq(x ,  t)  =  H(t  —  \x  —  x0|)/2, 
where  H(z )  is  the  Heaviside  function,  i.e., 


433 

434 

435 

436 

437 

438  AQ5 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 


H(z)  = 


(  0,  z  <  0 
\  1,Z>  0. 


Hence,  u(0,t)  =  g(t)  =  H(t  —  l)/2  for  1  <  t  <  1.4;  and  at  451 
t  =  1.4,  the  return  wave  from  the  block  hits  the  observation  452 
point  {x  =  0}  for  the  first  time.  453 

The  measured  data  are  also  challenging  to  process  since  454 
they  arise  from  oblique  illumination,  and  the  exact  location  455 
and  the  amplitudes  of  the  incident  pulses  were  not  known.  In  456 
addition,  a  comparison  of  Fig.  4(b)  with  Fig.  5(b),  (d),  and  (f)  457 
shows  that  the  measured  data  are  highly  oscillatory,  which  are  458 
unlike  their  simulated  counterparts.  Consequently,  we  applied  459 
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Fig.  4.  (a)  Function  er  (target)  =  4;  note  that  er  (background)  =  1.  (b)  u(0,  t)  =  g(t)  for  0  <  t  <  3.0.  The  source  is  located  at  xq  =  — 1,  and  the  first 

backscatter  return  is  therefore  shown  at  approximately  t  =  2.4  with  “ringing”  determined  by  interference  of  multiply  scattered  waves  between  the  two  boundaries 
of  the  block.  Computations  were  performed  using  the  software  package  WavES  [24]. 
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Fig.  5.  Three  targets  and  their  associated  measured  data.  The  ground  is  dry  sand  with  3  <  er  <  5  [21],  [22].  The  information  shown  in  (a),  (c),  and  (e)  were 
only  provided  after  computations  were  made,  (a)  Depicts  a  bush  that  was  located  on  a  road,  which  generated  background  clutter,  (b)  Scaled  experimental  data  for 
(a),  where  the  horizontal  axis  represents  time  in  nanoseconds  having  a  time  step  of  0.133  ns;  and  the  vertical  axis  is  the  amplitude  of  the  measured  voltage  at  the 
detector,  (c)  Wooden  stake,  (d)  Scaled  experimental  data  for  (c).  (e)  Metal  box  buried  in  dry  sand,  (f)  Scaled  experimental  data  for  (e).  The  mismatch  between 
experimental  and  simulated  data  [see  Fig.  4(b)]  is  evident. 
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Fig.  6.  Superimposed  preprocessed  data  for  all  five  cases  under  consideration. 
The  upward-looking  peak  corresponds  to  the  plastic  cylinder  (see  Table  I). 

460  an  intuitively  reasonable  data  preprocessing  procedure,  which 

461  remained  totally  unbiased  since  it  was  applied  to  blind  data  sets. 

462  The  idea  of  this  procedure  is  to  make  the  data  more  similar  to 

463  that  shown  in  Fig.  4(b).  Previously,  a  similar  procedure  was 

464  reported  for  transmitted  data  in  [6],  [7],  and  [12].  We  have 

465  considered  two  cases. 

466  Case  1.  Suppose  that  the  target  is  located  above  the  ground. 

467  In  this  case 

er (target)  >  sr (background)  =  £r(air)  =  1.  (3.3) 

468  Fig.  4(a)  and  (b)  shows  that,  in  this  case,  the  backscattering 

469  signal  should  be  basically  one  downward-looking  peak. 

470  Therefore,  we  have  selected  on  the  experimental  curve  the 

47 1  first  downward-looking  peak  with  the  largest  amplitude.  As 

472  to  the  rest  of  the  experimental  curve,  it  was  set  to  zero. 

473  Hence,  we  work  only  with  the  selected  peak. 

474  Case  2.  Suppose  that  the  target  is  buried  in  the  ground.  In  this 

475  case,  we  cannot  claim  the  validity  of  (3.3).  On  the  other 

476  hand,  our  numerical  simulations  (not  shown  here)  have 

477  demonstrated  that,  if  er (target)  <  er (background),  then 

478  in  the  analog  of  Fig.  4(b),  the  peak  would  look  upward. 

479  Therefore,  in  this  case,  we  have  selected  on  the  experimen- 

480  tal  curve  of  the  first  peak  with  the  largest  amplitude  to  work 

481  with  initially. 

482  We  were  provided  with  five  data  sets.  Fig.  6  shows  superim- 

483  posed  preprocessed  curves  for  all  five  targets  we  have  worked 

484  with.  The  only  peak  that  looks  upward  is  the  one  for  the  plastic 

485  cylinder  buried  in  soil  since  its  dielectric  constant  was  less 

486  than  that  of  the  soil  (see  Fig.  6).  We  stress  once  again  that 

487  nothing  was  known  in  advance  about  the  dielectric  constants 

488  of  targets.  Therefore,  the  choice  of  the  upward-looking  peak 

489  in  the  case  of  the  plastic  cylinder  was  unbiased  and  was  done 

490  only  using  the  aforementioned  rule.  The  measured  amplitude 

491  for  each  case  was  on  the  order  of  105.  This  is  well  above  the 

492  amplitude  in  Fig.  4(b).  Thus,  all  signals  were  preprocessed 

493  first  (as  previously  described)  and  multiplied  by  the  scaling 


Fig.  7.  Graphs  of  the  function  m(0,  s)  =  w( 0,  s)  —  mo(0,  s)  for  s  E  [1,  5] 
for  the  LT  of  the  computationally  simulated  data  in  Fig.  4(b)  and  of  the 
preprocessed  signal  for  the  bush  (see  Fig.  6).  The  signal  in  Fig.  6  for  the  bush 
was  multiplied  by  10-7.  Minimal  and  maximal  values  of  the  function  w(0:s) 
are  similar  for  both  curves.  A  similar  observation  was  made  for  four  other 
targets  we  have  worked  with. 

number  SN  =  1CT7  next.  Consider  the  LT  of  the  simulated  494 
data  shown  in  Fig.  4(b)  and  then  the  pre-processed  signal  for  the  495 
bush  (see  Fig.  6)  and  multiply  it  by  1CT7.  Fig.  7  depicts  super-  496 
imposed  graphs  of  the  function  w( 0,  s )  =  rc(0,  s )  —  u>o(0,  s )  497 
for  s  E  [1,  5]  for  both  cases.  One  can  see  that  maximum  and  498 
minimum  values  of  both  curves  are  approximately  the  same.  499 
We  initially  used  SN  =  1(T6,  SN  =  1(T7,  and  SN  =  1(T8.  500 
Only  for  SN  =  10-7,  the  maximum  and  minimum  values  of  501 
functions  w( 0,  s )  for  sG  [1,5]  of  both  aforementioned  curves,  502 
i.e.,  the  one  for  the  LT  of  the  function  depicted  in  Fig.  6  (bush),  503 
being  multiplied  by  10-7,  and  the  one  for  the  LT  of  the  function  504 
in  Fig.  4(b),  were  approximately  the  same.  On  the  other  hand,  505 
those  minimal  and  maximal  values  were  quite  different  from  the  506 
values  of  the  LT  of  the  function  in  Fig.  4(b)  for  SN  =  10-6  and  507 
SN  =  1CT8.  Using  SN  =  10-7,  which  is  based  on  the  data  for  508 
the  bush,  we  have  multiplied  the  other  four  preprocessed  signals  509 
(see  Fig.  6)  by  10-7  and  observed  a  similar  behavior  for  the  four  510 
other  targets.  For  the  case  of  the  inverted  peak  in  Fig.  6,  we  511 
compared  \w(0,  s)|  for  it  with  w( 0,  s )  for  the  aforementioned  512 
simulated  data.  Note  that  the  signals  shown  in  Fig.  6  are  not  yet  513 
multiplied  by  the  scaling  number.  After  multiplying  these  data  514 
by  the  scaling  factor  10-7,  then  for  each  set  of  experimental  515 
data,  we  took  the  resulting  curve  as  the  function  u(0,t)  —  516 
^o(0,t)  :=  g(t)  —uo(0,t).  Next,  we  worked  only  with  this  517 
function  as  the  data,  using  the  aforementioned  algorithm.  For  518 
simple  isolated  targets,  these  steps  of  data  preprocessing  are  519 
justified,  given  the  accuracy  of  the  results  obtained  upon  a  520 
posteriori  inspection.  For  more  complex  target  volumes,  a  more  521 
sophisticated  analysis  of  sets  of  time  histories  will  be  necessary.  522 
The  data  sets  were  processed,  and  the  targets  are  illustrated  in  523 
Fig.  5.  If  we  compare  the  highly  oscillatory  curves  of  Fig.  5(b),  524 
(d)  and  (f),  one  can  see  that  these  backscatter  time  histories  or  525 
signatures  are  qualitatively  quite  similar  in  appearance.  Their  526 
oscillatory  nature  is  due  to  the  specific  carrier  frequency  and  527 
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TABLE  I 

Computed  Values  for  R,  the  Relative  Dielectric  Constant  in  (3.1),  Based  on  Blind  Processing  of  Measured 
Backscatter  Data  From  Five  Different  Targets.  Here,  A  Means  Air,  and  B  Means  Dry  Sand 


Target 

A/B 

R 

er  (backgr) 

er  (target) ,  calc. 

£r  (target) ,  published. 

Figure  3.3-(a) 

n/a 

3.8 

i 

3.8 

4  (known) 

Bush 

A 

6.5 

i 

6.5 

3  to  20  [10] 

Wood  stake 

A 

3.8 

i 

3.8 

2  to  6  [21] 

Metal  box 

B 

3.8 

3  to  5  [21] 

11.4  to  19 

10  to  30  (3.2) 

Metal  cylinder 

B 

4.3 

3  to  5  [21] 

12.9  to  21.4 

10  to  30  (3.2) 

Plastic  cylinder 

B 

0.4 

3  to  5  [21] 

1.2  to  2 

1.1  to  3.2  [21,  22] 

528  finite  bandwidth  of  the  pulsed  radiation,  whereas  the  simulated 

529  data  assume  an  idealized  pulse.  For  these  simple  targets,  we 

530  allow  the  aforementioned  preprocessing  step  to  force  a  cor- 

531  respondence  between  the  two  in  order  to  identify  the  earliest 

532  return  from  the  boundary  of  the  target  and  determine  its  relative 

533  amplitude.  Based  on  this,  the  inversion  algorithm  can  determine 

534  a  reliable  estimate  of  that  target’s  actual  permittivity.  In  addi- 

535  tion,  we  have  conducted  a  limited  sensitivity  study  with  respect 

536  to  the  scaling  factor.  Specifically,  we  took  SN  =  0.8  •  1CT7 

537  and  SN  =  1.2  •  10“ 7  for  all  five  targets,  which  are  variations 

538  of  20%  of  the  scaling  number.  In  five  out  of  five  cases  of 
AQ6  539  experimental  data,  we  have  worked  with  values  of  R  kept  within 

540  tabulated  limits  (see  Table  I)  when  these  variations  of  SN 

541  were  tried.  An  optimal  value  of  SN  might  be  determined  via 

542  a  comparison  of  values  of  R  :=  R(SN)  with  measured  values 

543  for  a  few  known  targets.  At  present,  we  have  concentrated  on 

544  reconstructing  a  real  parameter  that  describes  the  permittivity 

545  of  target  features;  and  metal  objects  have  been  images  simply 

546  having  a  very  large  relative  permittivity.  We  note  that  there  is 

547  no  reason  why  a  conductivity  term  could  not  be  incorporated 

548  into  the  algorithm. 

549  In  addition  to  high  oscillations  of  the  data,  we  have  faced 

550  two  more  uncertainties.  First,  we  did  not  know  where  the 

551  time  t  =  0  is  on  our  data.  Second,  we  did  not  know  where 

552  the  actual  location  of  the  source  xo  is.  This  means  that  it  is 

553  impossible  to  determine  the  location  of  the  target.  Hence,  for 

554  computational  purposes,  we  have  arbitrarily  assigned  t  =  0  to 
AQ7  555  be  a  fixed  location  1  ns  off  to  the  left  from  the  beginning  of  the 

556  largest  amplitude  peak  and  xq  :=  —1,  knowing  that  we  have 

557  independent  GPS  data  to  better  fix  absolute  ranges  should  we 

558  need  that  information.  Our  primary  objective  here  is  to  confirm 

559  the  quantitative  accuracy  of  the  estimates  of  the  dielectric 

560  constant  of  each  of  the  targets,  i.e.,  to  accurately  image  the  ratio 

561  R  in  (3.1). 

562  The  derivative  of  the  LT  of  the  preprocessed  data  was  corn- 

563  puted  for  0  <  s  <  12  with  a  step  size  of  As  =  0.05.  Since 

564  the  calculation  of  the  derivative  of  noisy  data  is  an  ill-posed 

565  problem,  we  have  used  the  following  well-known  formula  for 

566  the  calculation  of  the  derivative  of  the  LT: 

oo 

<p'(s)  -  dsw0(0,  s)  =  -  J  ( g(t )  -  uq (0,  t))  te~stdt.  (3.4) 
o 

567  Since  for  all  targets  the  function  g(t)  —  ^o(0,  t)  =  0  for  t  >  2 

568  (see  Fig.  6),  then  the  integration  in  (3.4)  is  actually  carried  for 

569  0  <  t  <  2.  We  then  define  boundary  conditions  for  functions 

570  qn,k  f°r  each  n,  and  R  is  calculated  by  the  aforementioned 

571  algorithm. 


In  Fig.  8(a)  and  (f),  we  regard  R  as  the  maximal  amplitude  of  572 
the  calculated  peak.  We  first  verified  that  the  algorithm  provides  573 
a  good  estimate  for  R  using  simulated  data.  For  the  block  in  574 
Fig.  4(a),  we  obtain  the  1-D  image  shown  in  Fig.  8(a),  which  575  aqs 
was  found  to  be  er  =  3.8,  which  is  very  close  to  the  known  576 
value  of  4.  Next,  we  have  calculated  images  from  experimental  577 
data.  In  addition  to  Fig.  5(a),  (c),  and  (e),  we  have  also  imaged  578 
two  more  cases,  namely,  a  plastic  cylinder  and  a  metal  cylinder,  579 
which  are  both  buried  in  the  ground  with  schematics  similar  580 
with  the  one  in  Fig.  5(e).  Fig.  8(b)— (f)  displays  our  calculated  581 
images  for  all  five  targets.  582 

Dielectric  constants  were  not  measured  when  the  data  were  583 
collected.  Therefore,  we  have  compared  computed  values  of  584 
dielectric  constants  with  those  listed  in  tables  [21],  [22].  Note  585 
that  these  tables  often  provide  a  range  of  values  rather  than  586 
exact  numbers;  but  given  this  caveat,  the  calculated  results  587 
for  these  materials  are  well  within  the  range  of  expectations  588 
(see  Table  I).  589 

IV.  Conclusion  590 

We  have  described  a  new  method  for  recovering  quanti-  591 
tatively  reliable  estimates  of  target’s  material  properties  (di-  592 
electric  constants)  from  backscatter  field  measurements.  The  593 
method  is  an  inverse  scattering  algorithm  based  on  a  rigorously  594 
formulated  CIP.  The  numerical  method  is  constructed  to  ensure  595 
global  convergence,  and  therefore,  it  avoids  stagnation  at  erro-  596 
neous  solutions  for  images  of  target  permittivity  distributions.  597 
Furthermore,  the  method  requires  no  prior  knowledge  of  the  598 
inhomogeneities  present  in  the  target  volume.  These  properties  599 
are  rigorously  guaranteed.  The  authors  are  unaware  of  alterna-  600 
tive  numerical  methods  with  similar  characteristics  for  the  case  601 
of  the  CIPs  making  use  of  such  limited  data.  602 

The  approach  was  evaluated  here  using  data  provided  by  the  603 
ARL  from  a  forward-looking  radar  system  without  any  prior  604 
knowledge  of  the  targets  being  used.  The  data  were  measured  605 
using  oblique  incidence  and  with  unknown  source  locations,  606 
and  thus,  some  assumptions  were  made  to  provide  the  necessary  607 
inputs  for  the  algorithm.  The  procedure  first  estimates  a  solution  608 
that  has  defined  error  given  the  quality  of  the  data  but  which  609 
is  guaranteed  to  be  reliable.  To  simplify  matters,  only  images  610 
of  dielectric  constants  were  recovered  in  order  to  validate  the  611 
quantitative  accuracy  of  the  approach.  Data  sets  were  prepro-  612 
cessed,  and  a  downrange  permittivity  profile  was  calculated.  613 
If  the  angular  spread  of  backscatter  time  histories  would  be  614 
measured,  then  its  additional  processing  would  provide  a  3-D  615 
image  with  a  high  spatial  resolution,  despite  the  use  here  of  a  616 
single  source  point  (see  [7,  Fig.  6.3]).  617 
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Fig.  8.  Calculated  images  of  targets.  The  ratio  R  in  (3.1)  is  regarded  as  the  maximal  amplitude  of  the  imaged  peak,  (a)  Image  for  computationally  simulated  data 
as  a  verification  of  the  accuracy  of  our  algorithm.  The  rectangular  block  and  the  curve  are  true  and  computed  profiles  of  the  dielectric  constant,  respectively.  The 
computed  target/background  contrast  R  =  3.8,  which  corresponds  to  a  5%  of  error,  (b)  Image  of  the  bush  [see  Fig.  2(a)].  The  calculated  er  (bush)  =  6.5,  which 
is  in  the  range  of  tabulated  values  3  <  er  <  20  [10].  (c)  Image  of  the  wood  stake  [see  Fig.  4(c)].  The  calculated  er(wood  stake)  =  3.8  [10].  (d)  Image  of  the 
buried  metal  box  [see  Fig.  5(e)].  The  calculated  R  =  3.8.  Since  the  background  was  dry  sand  with  3  <  er  (dry  sand)  <  5  [21],  then  the  calculated  er  (metal  box) 
is  between  11.4  and  19.  This  is  within  the  range  [see  (3.2)]  of  appearing  dielectric  constants  of  metallic  targets,  (e)  Calculated  image  of  the  buried  metal  cylinder. 
The  calculated  ratio  R  =  4.3.  Similarly  with  (d),  we  conclude  that  the  calculated  value  of  er  (metal  cylinder)  is  between  12.9  and  21.4.  This  is  again  within  the 
range  [see  (3.2)]  of  appearing  dielectric  constants  of  metallic  targets,  (f)  Calculated  image  of  the  buried  plastic  cylinder.  The  calculated  ratio  R  =  0.4.  Similarly 
with  (d),  we  conclude  that  the  calculated  value  of  the  dielectric  constant  er  (plastic  cylinder)  is  between  1.2  and  2.5,  which  is  again  within  the  range  of  tabulated 
values  for  plastic  [21],  [22]. 
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618  Since  the  dielectric  constants  of  the  targets  were  not  actually 

619  measured  in  the  ARL  experiments,  then  the  best  one  can  do 

620  is  compare  retrieved  parameters  with  tabulated  values.  Table  I 

621  shows  the  computed  relative  permittivities  of  targets.  It  is 

622  clearly  shown  that  all  five  targets  fall  well  within  expected 

623  tabulated  limits  for  the  materials  in  question,  despite  the  fact 

624  that  no  prior  knowledge  whatsoever  was  employed.  We  further 

625  emphasize  that  these  results  were  obtained  despite  a  very  lim- 

626  ited  information  content,  large  noise  in  the  data,  and  significant 

627  discrepancies  between  experimental  and  simulated  data.  We  can 

628  therefore  conclude  that  these  results  point  toward  the  validity  of 

629  our  mathematical  model.  The  fact  that  regardless  of  limitations 

630  of  the  method,  we  consistently  got  results,  which  only  later 

631  were  found  to  fall  well  within  tabulated  limits,  points  toward 

632  a  great  degree  of  robustness  of  this  algorithm. 

633  The  purpose  of  estimating  the  dielectric  constant  is  to  provide 

634  one  extra  piece  of  information  about  the  target.  Up  to  this  point, 

635  most  of  the  radar  community  has  solely  relied  on  the  intensity 

636  of  the  radar  image  for  doing  detection  and  discrimination.  It  is 

637  anticipated  that,  when  the  intensity  information  is  coupled  with 

638  the  new  dielectric  information  this  method  provides,  algorithms 

639  can  be  then  designed  that  will  provide  better  performance  in 

640  terms  of  probability  of  detection  and  false  alarm  rates.  Finally, 

641  we  repeat  that  the  results  presented  in  this  paper  are  primarily 

642  being  used  as  a  vehicle  to  illustrate  this  powerful  inverse 

643  scattering  algorithm  method  and  its  ability  to  recover  dielectric 

644  properties  of  targets  from  experimental  data  collected  by  the 

645  forward-looking  radar  of  the  ARL.  Detailed  studies  making 

646  use  of  larger  experimental  data  sets  from  more  complex  3-D 

647  scattering  objects  are  necessary,  and  the  authors  will  report  on 

648  this  in  the  near  future. 
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1  Quantitative  Image  Recovery  From  Measured  Blind 

2  Backscattered  Data  Using  a  Globally 

3  Convergent  Inverse  Method 

4  Andrey  V.  Kuzhuget,  Larisa  Beilina,  Michael  V.  Klibanov,  Anders  Sullivan,  Lam  Nguyen,  and  Michael  A.  Fiddy 


5  Abstract — The  goal  of  this  paper  is  to  introduce  the  application 

6  of  a  globally  convergent  inverse  scattering  algorithm  to  estimate 

7  dielectric  constants  of  targets  using  time-resolved  backscattering 

8  data  collected  by  a  U.S.  Army  Research  Laboratory  for  war  d- 

9  looking  radar.  The  processing  of  the  data  was  conducted  blind,  i.e., 

10  without  any  prior  knowledge  of  the  targets.  The  problem  is  solved 

11  by  formulating  the  scattering  problem  as  a  coefficient  inverse 

12  problem  for  a  hyperbolic  partial  differential  equation.  The  main 

13  new  feature  of  this  algorithm  is  its  rigorously  established  global 

14  convergence  property.  Calculated  values  of  dielectric  constants  are 

15  in  a  good  agreement  with  material  properties,  which  were  revealed 

16  a  posteriori . 

17  Index  Terms — Experimental  data,  inverse  scattering,  quantita- 

18  tive  imaging,  remote  sensing. 

19  I.  Introduction 

20  A  FUNDAMENTAL  problem  in  remote  sensing  is  the 

21  processing  of  scattered  field  data  from  strongly  scattering 

22  penetrable  targets.  Multiple  scattering  renders  this  problem  ex- 

23  tremely  difficult  to  solve,  it  being  ill  conditioned  with  additional 

24  questions  of  uniqueness  and,  the  most  difficult,  nonlinearity  to 

25  contend  with.  In  practice,  limited  noisy  data  typically  require 

26  that  some  physical  models  be  assumed,  from  which  one  hopes 

27  to  extract  meaningful  and  preferably  quantitative  information 

28  about  the  target  in  question.  A  number  of  recent  publications 

29  by  Beilina  and  Klibanov  [3]— [8]  and  by  Klibanov  et  al.  [12], 

30  [14]— [16]  have  led  to  a  new  approach  to  address  this  important 

31  topic.  This  numerical  method  was  originally  developed  for 

32  some  multidimensional  coefficient  inverse  problems  (MCIPs) 

33  for  a  hyperbolic  partial  differential  equation  (PDE)  using  data 
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from  only  a  single  location  of  either  a  point  source  or  from  34 
a  single  direction  of  an  incident  plane  wave.  In  particular,  in  35 
[14],  that  method  was  extended  from  the  3-D  case  to  the  1-D  36 
case.  Thus,  that  1-D  version  of  [14]  is  used  here  to  work  with  37 
the  experimental  data.  The  illuminating  field  is  pulsed  in  time,  38 
and  the  time  history  of  the  backscattering  from  the  illuminated  39 
target  volume  constitutes  the  measured  data  that  are  processed  40 
by  this  algorithm.  The  authors  are  unaware  of  other  groups  41 
working  on  MCIPs  using  data  acquired  from  a  single  source  42 
location.  However,  the  single  measurement  case  is  clearly  the  43 
most  practical  one,  particularly  for  military  applications.  In-  44 
deed,  because  of  many  dangers  on  the  battlefield,  the  number  45 
of  measurements  should  be  minimized.  46 

The  algorithm  in  the  aforementioned  cited  publications  com-  47 
putes  values  for  the  spatial  distribution  of  the  dielectric  con-  48 
stants  of  objects  within  the  target  volume.  It  is  important  to  49 
stress  that  this  algorithm  requires  neither  no  prior  knowledge  50 
of  what  might  exist  in  the  target  volume  nor  a  prior  knowledge  51 
of  a  good  first  guess  about  the  solution.  There  is  a  rigorous  guar-  52 
antee  that  this  algorithm  globally  converges  (see  mathematical  53 
details  in  [7],  [14],  [16],  and  [17]).  Because  of  the  global  con-  54 
vergence  property,  estimates  of  spatially  distributed  dielectric  55 
constants  are  reliable  and  systematically  improve  with  more  56 
measured  and  less  noisy  data.  The  theory  of  the  aforementioned  57 
cited  publications  rigorously  guarantees  that  this  numerical  58 
method  delivers  a  good  approximation  to  the  exact  solution  59 
of  an  MCIP  without  any  a  priori  information  about  a  small  60 
neighborhood  of  the  exact  solution  as  long  as  iterations  start  61 
from  the  so-called  “first  tail  function”  V${x),  which  can  be  62 
easily  computed  using  available  boundary  measurements  (see  63 
(2.27)-(2.29)  in  Section  II-C).  In  addition,  it  is  in  this  sense  64 
that  we  use  the  term  “global  convergence"  of  the  algorithm.  65 
The  common  perception  of  the  term  "global  convergence"  is  66 
that  one  can  start  from  any  point  and  still  get  the  solution,  but  67 
we  stress  that  we  actually  start  not  from  any  point  but  rather  68 
from  the  function  Vo(x ),  which  can  be  easily  computed  from  69 
the  boundary  data  (see  (2.27)-(2.29)  in  Section  II-C).  70 

It  is  well  known  that  least  squares  functionals  for  MCIPs  71 
suffer  from  multiple  local  minima  and  ravines.  Hence,  local  72 
convergence  of  numerical  methods  to  incorrect  estimates  will  73 
occur  unless  an  initial  guess  that  is  close  to  the  true  solution  is  74 
used.  Such  a  guess  is  rarely  available  in  most  applications.  In  75 
contrast,  our  algorithm  does  not  use  a  least  squares  functional,  76 
and  hence,  it  is  free  from  the  problem  of  local  minima.  Instead,  77 
this  algorithm  relies  on  the  structure  of  the  differential  operator  78 
of  the  wave-like  PDE.  79 
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80  Prior  to  the  work  reported  here,  a  major  focus  by  the 

81  U.S.  Army  Research  Laboratory  (ARL)  had  been  on  the  de- 

82  velopment  of  image  processing  techniques  [19]  that  would 

83  improve  radar  images,  which  is  through  postprocessing  tech- 

84  niques  rather  than  through  the  application  of  inverse  scattering 

85  methods.  By  incorporating  more  physics  of  the  target- wave 

86  electromagnetic  response  into  the  data  processing,  one  can 

87  greatly  improve  target  detection  and  identification.  Present  data 

88  processing  provides  an  electromagnetic  field  brightness  or  an 

89  intensity  map  of  the  target  volume,  which  need  not  relate 

90  in  a  simple  fashion  to  the  scattering  structures  themselves. 

91  Our  method  estimates  dielectric  constants  of  targets,  which 

92  obviously  adds  an  important  new  dimension  to  the  interpre- 

93  tation  of  data  acquired  by  the  radar  system  since  this  allows 

94  specific  bounds  on  the  dielectric  properties  of  a  feature  in 

95  the  target  volume,  which  can  help  identify  its  likely  material 

96  properties.  Since  no  prior  knowledge  is  required,  the  measured 

97  data  were  processed  by  Kuzhuget,  Beilina,  Klibanov,  and  Fiddy 

98  in  the  most  challenging  scenario,  i.e.,  without  any  knowledge 

99  of  the  actual  target  structures  and  their  dielectric  properties. 

100  Once  this  had  been  done,  Sullivan  and  Nguyen  compared  a 

101  posteriori  the  image  estimates  with  the  actually  known  material 

102  characteristics. 

103  We  draw  attention  to  the  fact  that  this  algorithm  has  been 

104  used  with  forward-scattered  data  from  experiments.  These 

105  results  were  previously  reported,  which  are  also  in  a  blind 

106  experiment  (see  [12,  Tables  5  and  6]  and  [7,  Tables  5.5  and 

107  5.6]).  In  this  case,  the  images  in  [12]  were  further  improved 

108  and  presented  in  a  follow-up  publication  [6]  using  the  adaptivity 

109  technique  of  [1],  [2],  [4],  [5],  and  [7]. 

110  In  Section  II,  we  outline  the  basic  steps  in  the  underlying 

111  theory  upon  which  the  new  algorithm  is  based.  In  Section  III, 

112  we  formulate  the  global  convergence  theorem.  In  Section  IV, 

113  we  outline  results  obtained  using  time-resolved  backscatter 

114  electric  field  measurements  collected  in  the  field.  Measure- 

115  ments  were  carried  out  by  a  forward-looking  radar  system  built 

116  and  operated  by  the  ARL.  The  data  were  noisy  and  limited,  and 

117  the  target  volumes  included  miscellaneous  sources  of  clutter. 

118  The  purpose  of  this  particular  radar  system  is  to  detect  and 

119  possibly  identify  shallow  explosive-like  targets. 


120 


II.  Theoretical  Background 


We  assume  that  the  constitutive  parameter  of  interest,  i.e.,  135 
mapping  the  target  volume,  is  a  relative  permittivity  er(x).  In  136 
other  words,  we  ignore  magnetic  effects  in  this  paper.  We  also  137 
assume  for  convenience  that  er(x)  =  1  outside  of  the  target  138 
volume,  which  is  x  £  (0, 1)  in  our  case.  We  assume  that  the  139 
source  xo  <  0  lies  outside  of  the  target  volume.  We  can  write  140 
the  forward  scattering  problem  as  141 

Er{pc)'Ujtt  — xx •>  X  £  M  (2.1) 

?i(x,0)=0,  ut(x,  0)  =  S(x  —  xo).  (2.2) 

The  subscripts  in  (2.1)  indicate  the  number  of  partial  derivatives  142 
with  respect  to  the  variable  indicated.  The  coefficient  inverse  143 
problem  (CIP)  is  to  recover  er(x ),  assuming  that  the  initial  144 
illuminating  pulse  is  known  and  that  we  measure  the  function  145 
g(t ),  i.e.,  146 


u(0,t)  =  g(t) 


(2.3) 


121  A.  Integral  Differential  Equation 

122  Since  we  were  given  only  one  time-resolved  experimental 

123  curve  per  target,  we  had  no  choice  but  to  use  a  1-D  mathemati- 

124  cal  model,  although  the  reality  is  3-D  (see  Section  III  for  some 

125  details  about  the  data  collection).  In  addition,  since  only  one 

126  component  of  the  electric  wave  field  was  both  transmitted  and 

127  measured,  we  model  the  scattering  process  with  one  wave-like 

128  PDE  rather  than  using  complete  Maxwell  equations.  We  stress 

129  that  the  method  is  designed  for  use  with  3-D  problems,  and 

130  one  would  normally  collect  data  with  co  polarization  and  cross 

131  polarization  in  order  to  capture  all  of  the  pertinent  information 

132  about  the  target.  Here,  we  simply  wish  to  show  the  steps 

133  employed  by  the  method  and  demonstrate  their  quantitative 

134  reconstruction  accuracy  given  noisy  measured  data. 


for  sufficiently  large  times  t  that  all  multiple  scattering  events  147 
within  the  target  volume,  which  can  produce  a  measurable  148 
signal  at  the  detector,  do  so.  Practically,  we  gate  the  radiation  149 
source  in  time;  and  since  the  Laplace  transform  (LT),  i.e.,  150 
w(x,s),  is  used  to  solve  this  CIP,  the  decay  e~st,s  >  0  of  151 
the  LT  kernel  further  limits  the  duration  of  the  measured  time  152 
history.  It  is  worth  pointing  out  that,  more  typically,  scattering  153 
data  would  be  measured  at  different  scattering  angles  for  fixed  154 
frequency  illumination  at  various  incident  angles.  One  can  155 
easily  appreciate  that  this  leads  to  the  acquisition  of  Fourier  156 
information  about  the  target  or  the  secondary  source  function,  157 
depending  upon  the  extent  of  the  multiple  scattering;  and  once  158 
one  has  sufficient  data,  a  reasonable  estimate  of  the  target  159 
properties  becomes  possible.  By  taking  measurements  in  the  160 
time  domain,  one  can  see  that  this  is  essentially  simultane-  161 
ously  gathering  information  in  a  transform  space  from  many  162 
illumination  frequencies.  The  Laplace  and  Fourier  transforms  163 
provide  complimentary  representations  of  the  target  in  terms  of  164 
moments  or  modes,  respectively.  165 

The  LT  is  166 

03 

M*,s)  =  ju(x,t}e-<,t:=Cu,  »>s=conS,.>0  (2.4) 

0 

and  we  assume  that  the  so-called  pseudofrequency  s  >  167 
s(er(x))  :=  s  is  sufficiently  large.  This  gives  [7]  168 

wxx  ~  s2er(x)w  =  —  S(x  —  xq),  x  £  M  (2.5) 


lim  w(x,  s )  =  0. 

|^-0O3 


(2.6) 


Let 


169 


w(  0,  s)  =  ip(s)  =  Cg 


(2.7) 


be  the  LT  of  the  measured  function  g(t)  in  (2.3).  Since  er(x)  =  170 
1  for  x  <  0,  then,  using  (2.5)  and  (2.6),  one  can  prove  that,  in  171 
addition  to  the  function  rc(0,  s)  in  (2.7),  the  function  wx( 0,  s)  172 
is  also  known  as  (see  [17])  173 


Wx( 0,  s)  =  Sip(s)  -  exp(sxo). 


(2.8) 
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174  Let  wo  (x,  s)  be  the  solution  of  the  problem  in  (2.5)  and  (2.6) 

175  for  the  case  of  the  uniform  background  er(x)  =  1.  Then 

exp  (— s\x  —  #o|) 


WQ(X,S)  = 


2s 


r(x,s)  =  4ln  (  —  (x,s)  )  . 

Sz  \Wo 


qxx  +  2s2qxrx  +  2srl~2sqx-2rx=0, 


V (x,  s )  :=  V (x)  =  r(x,  s)  =  4  ln  ( 

s  \Wo(x,S) 


Qxx  2s  qx 

—  2 sqx  +  2 


ir  +  2s 


J  qx(x,r)d 

S 

s 

J  qx(x,T)dT 

S 

S 

+  2s2qxVx  -  4sVx  J  qx(x,  r)d 

+  2s{Vxf-2Vx  =0, 
x  G  (0, 1);  s  G  [s,  s] 
q(0,s)  =  ipo(s),  qx(0,s)  =  i(>1(s) 
qx( l,s)sO,  s  G  [s,  s] 


J  qX{x,T)d7 


(2.9) 


176  When  implementing  the  algorithm,  given  the  assumption  of  a 

177  uniform  normalized  er(x)  =  1  outside  of  the  target  volume,  we 

178  consider  the  function 


(2.10) 


179  Since  the  source  xo  <  0,  then  the  function  r(x,  s )  is  the  solution 

180  of  the  following  equation  in  the  interval  (0,  1): 

rXx  +  s2r2x  -  2srx  =  er(x)  -  1,  iG(0,l).  (2.11) 

181  In  addition,  by  (2.7)  and  (2.8) 

r(0,s)  =ip0(s),  rx(0,  s)  =  ipi(s)  (2.12) 

ln</?(s)  -  ln(2s)  x0 

lPo{s)  ~ - ~2 - ^  GT 

sz  s 

2  esx° 

<Pi(*)= - 2-TT-  (2'13) 

S  s2ip(s) 

182  The  idea  now  is  to  eliminate  the  unknown  coefficient  er(x) 

183  from  (2.11)  via  differentiation  with  respect  to  pseudofre- 

184  quency  s.  Differentiating  (2.11)  with  respect  to  s  and  denoting 

185  q(x,  s )  =  dsr(x ,  s),  we  obtain 


x  G  (0, 1).  (2.14) 


186  We  now  need  to  express  in  (2. 14)  the  function  r  via  the  function 

187  q.  We  have 

s 

>■(*.»)  =  -  j  q(x,T)dT  +  V(x,H)  (2.15, 


188  where  V{x)  :=  V  (x,  s)  is  referred  to  as  the  tail  function,  which 

189  is  small  in  practice  for  large  positive  s.  Here,  the  truncation 

190  pseudofrequency  s  serves  as  a  regularization  parameter.  The 

191  exact  formula  for  V  (x)  is 


Fig.  1.  Flowchart  of  the  algorithm. 

where  functions  ^o( <5 )  =  <Po(s)  and^i(s)  =  ipi(s)  are  derived  194 
from  (2.13).  The  condition  qx(l,  s)  =  0  can  be  easily  derived  195 
from  (2.6)  since  er{pc)  =  1  outside  of  the  interval  (0,  1).  196 

In  (2.17)  and  (2.18),  both  functions  q(x,s)  and  V(x)  are  197 
unknown.  The  reason  why  we  can  approximate  both  of  them  198 
is  that  we  find  updates  for  q(x,s)  via  inner  iterations  exploring  199 
(2.17)  and  (2.18)  inside  of  the  interval  (0,  1).  At  the  same  time,  200 
we  update  the  tail  function  V  (x)  via  outer  iterations  exploring  201 
the  entire  real  line  M.  In  short,  given  an  approximation  for  V (x),  202 
the  algorithm  updates  q  and  then  updated  for  er(x).  Next,  the  203 
forward  problem  in  (2.5)  and  (2.6)  is  solved  for  the  function  204 
w(x,s)  for  s  =  s.  Next,  the  tail  function  V (x)  is  updated  using  205 
(2.16).  This  might  seem  reminiscent  of  the  steps  in  algorithms  206 
such  as  the  modified  gradient  inverse  scattering  technique  [20] ;  207 
but  we  emphasize  that,  unlike  our  case,  such  methods  have  no  208 
global  convergence  properties.  209 


B.  Iterative  Process 


210 


(2.16) 


192  Substituting  (2.15)  in  (2.14),  we  obtain  the  following  nonlinear 

193  integral  differential  equation: 

-i  2 


We  now  outline  the  formulation  of  our  algorithm  and  the  211 
iterative  process  (see  details  in  [7],  [14],  [16],  and  [17];  see  212 
Fig.  1).  Unlike  computationally  simulated  data  in  [14],  we  213  AQ2 
do  not  use  prior  knowledge  of  the  function  g(l,s)  on  the  214 
transmitted  edge  since  this  function  is  unknown  to  us.  We  have  215 
observed  in  our  computational  experiments  that  the  knowledge  216 
of  q(l,s)  only  affects  the  accuracy  of  the  calculation  of  the  217 
location  of  the  target,  but  it  does  not  affect  the  accuracy  of  the  218 
computed  target/background  contrast.  Here,  we  are  interested  219 
only  in  that  contrast  (see  Section  III).  Since  er(x)  =  1  for  x  >  1  220 
and  xo  <  0,  then  one  can  easily  derive  from  equations  (2.5),  221 
(2.9),  and  (2.10)  that  dxq(  1,  s )  =  0.  222 

Consider  a  partition  of  the  interval  [s,  s]  into  N  small  subin-  223 
tervals  with  the  small  grid  step  size  h  >  0  and  assume  that  the  224 
function  q(x,  s)  is  piecewise  constant  with  respect  to  s.  Thus  225 


S  =  S tv  <  Sjv-1  <  •  •  •  <  So  =  5, 
q(x,s)  =qn(x),  for  s  G  (sn,  sn_i]. 


Si- 1  —  Si  =  h 

(2.19) 


(2.17) 

(2.18) 


For  each  subinterval  (sn,  sn-i]  we  obtain  a  differential  equation  226 
for  the  function  qn(x).  We  assign,  for  convenience  of  notations,  227 
qo  :=  0.  Following  the  aforementioned  idea  of  a  combination  of  228 
inner  and  outer  iterations,  we  perform  for  each  n  inner  iterations  229 


230  with  respect  to  the  tail  function.  This  way,  we  obtain  functions 

231  </„,/,  and  Vn,k.  The  equation  for  the  pair  14, fc)  is 


n— 1 


Qn,k  1  A\ ,n/z  ^  ]  Qj  AijnVn,k  2^2, n  I  Qn,k 


3= 0 
^  n—1 


n— 1 


n—1 


~A2,nh2  {'El]  +2h^j+2A2tnV^k\h^j 


vi=° 
\  2 


i=o 


i=o 


V'  , 

i  v  n^ki 

xe(0,i) 

O 

Jl^ 

t—i 

3 

II 

1 

sn—  1 
/» 

(2.21) 


232  Here,  A\,n  and  A2,n  are  certain  numbers,  whose  exact  expres- 

233  sions  are  given  in  [3]  and  [7]. 

234  The  choice  of  the  first  tail  function  Vo(x)  is  described  in 

235  Section  II-C.  Let  n  >  1.  Suppose  that,  for  j  =  0, . . .  n  —  1, 

236  functions  qj{x)  and  Vj(x)  are  already  constructed.  We  now 

237  need  to  construct  functions  qn^  and  Vnjk  for  k  =  1, . . . ,  m. 

238  We  set  Vn,i(x)  :=  Vn-\{x).  Next,  using  the  quasi-reversibility 

239  method  (QRM)  (see  Section  II-C),  we  approximately  solve 

240  (2.20)  for  k  —  1  with  overdetermined  boundary  conditions  in 

241  (2.21)  and  find  the  function  qn^.  Next,  we  find  the  approxima- 

242  tion  for  the  unknown  coefficient  er{x)  via  the  following 

243  two  formulas: 

n—l 

rn, lW  =  -  hqn, i  -  h^qj  +  14, l,  x  G  [0,1]  (2.22) 

J=0 

4n,1)(X)  =  i  +  r"^)  +  4  [4,i  60] 2 

-2snr'nA(x),  x  e  [0, 1].  (2.23) 

244  Next,  we  solve  the  forward  problem  in  (2.5)  and  (2.6)  with 

245  er(x)  :=  £rn,1)(^),  s  :=s  and  find  the  function  wn^(x,s) 

246  this  way.  After  this,  we  update  the  tail  via  the  formula  in  (2.16), 

247  in  which  w(x,  s)  :=  wn^(x,  s).  This  way,  we  obtain  a  new  tail 

248  Vn^{x).  Similarly,  we  continue  iterating  with  respect  to  tails  m 

249  times.  Next,  we  set 

qn(x)  ■■=  qn,m(x),  Vn(x)  :=  Vn,m(x),  e^\x :)  := 

250  replace  n  with  n  +  1  and  repeat  this  process.  We  continue  this 

251  process  until  [15] 


either 


£(n)  _  ^O-1) 


MOT) 


<  10 


-5 


or  ||VJa(gnjfe)||L2(0jl)  >  105  (2.24) 


252  where  the  functional  Ja{qn,k )  is  defined  in  Section  II-C.  Here, 

253  the  norm  in  the  space  L2(0, 1)  is  understood  in  the  discrete 

254  sense.  In  the  case  when  the  second  inequality  in  (2.24)  is 

255  satisfied,  we  stop  at  the  previous  iteration,  taking  e (x)  as 

256  our  solution.  If  neither  of  two  conditions  in  (2.24)  is  not  reached 

257  at  n  :=  N,  then  we  repeat  the  aforementioned  sweep  over  the 

258  interval  [s ,  s] ,  taking  the  pair  (q]y(x),  Vn(x))  as  the  new  pair 

259  (qo(x),  V& (a)).  Usually,  at  least  one  of  the  conditions  in  (2.24) 

260  is  reached  either  on  the  third  or  on  the  fourth  sweep,  and  the 
AQ3  261  process  stops  then. 
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C.  Computing  Functions  qn,k{x)  and  Vq(x)  262 

At  first  glance,  it  seems  that,  for  a  given  tail  function  Vn^  (x),  263 
the  function  qn,k(x)  can  be  computed  as  the  solution  of  a  264 
conventional  boundary  value  problem  for  the  ordinary  differ-  265 
ential  equation  in  (2.20)  with  any  two  out  of  three  boundary  266 
conditions  in  (2.21).  However,  attempts  to  do  so  led  to  poor  267 
quality  images  (see  [14,  Remark  3.1]).  At  the  same  time,  the  268 
QRM  has  resulted  in  accurate  solutions  both  in  [14]  and  in  Test  269 
1  for  synthetic  data  (see  succeeding  discussion).  The  QRM  is  270 
well  designed  to  compute  least  squares  solutions  of  PDEs  with  271 
overdetermined  boundary  conditions,  such  as,  e.g.,  the  problem  272 
in  (2.20)  and  (2.21).  We  refer  to  [18]  for  the  originating  work  273 
about  the  QRM  and  to  [7],  [9],  [13],  [15],  and  [16]  for  some  274 
follow-up  publications.  275 

Let  L(qn,k)(x)  and  Pn^{x)  be  left-  and  right-hand  sides  of  276 
(2.20),  respectively.  In  our  numerical  studies,  L(qn^)(x)  and  277 
PU:k{x)  are  written  in  the  form  of  finite  differences.  Let  a  G  278 
(0, 1)  be  the  regularization  parameter.  The  QRM  minimizes  the  279 
following  Tikhonov  regularization  functional:  280 

Jaijln^k)  =  ||  Ln,k  ( qn,k )  Pn,k  IIl2(0,1)  ^ll^n,fc  II  JT2(0,1)  (2-25) 

subject  to  boundary  conditions  in  (2.21).  Here,  again  norms  281 
in  L2(0, 1)  and  in  the  Sobolev  space  H2( 0, 1)  are  understood  282 
in  the  discrete  sense.  The  functional  Ja(qn,k )  in  (2.25)  is  283 
quadratic.  Using  this  fact  and  the  tool  of  Carleman  estimates,  it  284 
can  be  shown  that  Ja(qn,k )  has  a  unique  global  minimum  and  285 
no  local  minima  [14],  [15],  [17].  We  find  that  global  minimum  286 
via  the  conjugate  gradient  method,  minimizing  with  respect  to  287 
the  values  of  the  function  qn^  at  grid  points.  We  have  used  288 
100  grid  points  in  the  interval  (0,  1).  The  step  size  in  the  s-  289 
direction  was  h  =  0.5.  The  8-interval  was  [s,  s]  =  [3, 12].  For  290 
each  n  =  1, . . . ,  N,  we  take  functions  qn^  for  k  =  1, ...  m,  291 
and  we  typically  choose  m  =  10.  The  reason  for  the  choice  292 
of  m—  10  is  that  numerical  experience  has  shown  that,  for  293 
each  of  the  n ,  tails  stabilize  at  k  «  10.  As  to  the  regularization  294 
parameter  a ,  we  have  found,  when  testing  synthetic  data,  that  295 
ot  =  0.04  is  the  optimal  one,  and  we  take  it  in  our  computations.  296 

We  note  that  we  determined  the  regularization  parameter  297 
when  testing  simulated  data.  These  data  were  for  the  target  298 
depicted  in  Fig.  7(a),  for  which  we  varied  the  regularization  299 
parameter  between  0.03  and  0.05.  The  resulting  images  for  300 
these  data  showed  only  an  insignificant  change.  We  also  var-  301 
ied  the  regularization  parameter  between  0.03  and  0.05  for  302 
the  experimental  data.  Again,  we  only  observed  insignificant  303 
changes,  which  lead  us  to  select  the  average  value  of  0.04.  304 
Although  the  regularization  theory  states  that  the  regularization  305 
parameter  should  depend  on  the  noise  level  in  the  data  [23],  we  306 
do  not  actually  know  the  noise  level  for  our  data.  Further,  for  307 
nonlinear  problems  (as  we  have),  this  dependence  is  claimed  308 
by  regularization  theory  only  for  the  limiting  case  of  a  relatively  309 
small  level  of  noise,  which  is  not  our  case.  In  our  computations  310 
using  measured  data,  one  works  with  some  level  of  noise,  which  311 
is  not  likely  to  be  small  and  is  unknown.  Therefore,  in  practice,  312 
when  applying  this  algorithm  to  experimental  data,  we  were  313 
guided  by  results  from  simulations  to  choose  a  value  for  the  314 
regularization  parameter.  If  we  had  prior  knowledge  about  some  315 
objects  in  the  target  volume,  then  we  would  choose  the  optimal  316 
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317  regularization  parameter  for  that  object.  Because  we  processed 

318  the  data  without  any  prior  knowledge  whatsoever  about  the 

319  objects,  we  chose  the  regularization  parameter  based  on  the 

320  simulated  data  processing,  and  fortunately,  our  answers  for  five 

321  out  of  five  targets  were  well  within  tabulated  limits. 

322  We  now  describe  an  important  step  in  choosing  the  first 

323  tail  function  Vq(x).  To  choose  it,  we  consider  the  asymptotic 

324  behavior  of  the  function  V(x,s)  in  (2.16)  with  respect  to  the 

325  truncation  pseudofrequency  s  -A  oo.  This  behavior  is  [14],  [17] 


V(x,s)  =  V-^ 

s 


0IP 


oo. 


326  We  truncate  the  term  0(1/ s2),  which  is  somewhat  similar  with 

327  the  defining  of  geometrical  optics  as  a  high-frequency  approx- 

328  imation  of  the  solution  of  the  Helmholtz  equation.  Hence,  we 

329  take 


V(x,s) 


Po(x) 


330  Since  q  =  dsr  and  V (x,  s )  =  r(x,  s ),  then 

Po(x) 


q(x,s)  = 


s2 


(2.26) 


331  Hence,  setting  in  (2.17)  s  :=  s  and  using  (2.26),  we  obtain  the 

332  following  approximate  equation  for  the  function  po(x): 


d 2 

^»(*)  =  o, 


X  £  (0, 1). 


(2.27) 


V  o(-r)  := 


(2.29) 


IV’o(s)  -  $S(S) I  <  <7,  IV’i(s)  -  V’iOOI  < 


for  s  E  Is,  s } 


346  where  functions  'ipo(s)  and  Vh  (s)  depend  on  the  function  g(t)  in 

347  (2.3)  via  (2.7),  (2.13)  and  (2.18);  and  functions  'iPq(s)  and  'ipl(s) 

348  depend  on  the  noiseless  data  g*  ( t )  in  the  same  way.  Let  h  E 

349  (0, 1)  be  the  grid  step  size  in  the  5-direction  in  (2.19);  let  yj a  = 

350  cr  and  h  =  max(cr,  h).  Let  Q  be  the  total  number  of  functions 

(n  k) 

351  Sr  ’  computed  in  the  aforementioned  algorithm.  Then,  there 

352  exists  a  constant  D  =  D(xo,  d,  s)  >  1  such  that,  if  the  numbers 

353  a  and  h  are  so  small,  that 

^  <  JJ2Q+2  (2.30) 


Fig.  2.  (a)  Schematic  diagram  of  the  forward-looking  radar  system  illumi¬ 

nating  a  dielectric  target,  (b)  Typical  measured  time  history  of  the  backscatter 
held,  (c)  Composite  of  unprocessed  returns  highlighting  the  dielectric  target 
(indicated  by  the  red  triangle),  (d)  Downrange  cut  of  the  permittivity  profile, 
which  the  new  algorithm  will  generate. 


then  the  following  estimate  is  valid: 

r(n,/c)  _ 


354 


L2(  0,1) 


<  hu 


(2.31) 


_(n,k) 


333  Boundary  conditions  for  po(x)  can  be  easily  derived  from 

334  (2.18)  and  (2.26)  as 

Po(0)  =  ~s2Ms),  Po(0)  =  ~s2Ms),  PbO)  =  0-  (2.28) 

335  We  find  an  approximate  solution  po,appr(x)  of  the  problem  in 

336  (2.27)  and  (2.28)  via  the  QRM,  similarly  with  the  aforemen- 

337  tioned  equation.  Next,  we  set  for  the  first  tail  function,  i.e., 

P0,appr  (x) 


338  A  simplified  formal  statement  of  the  global  convergence 

339  theorem  is  as  follows  (see  [7,  Th.  6.1]  for  more  details  and 

340  [7,  Th.  6.7]  for  the  3-D  case). 

341  Theorem  1:  Let  the  function  e*r(x)  be  the  exact  solution  of 

342  our  CIP  for  the  noiseless  data  g*(t)  in  (2.3).  Fix  the  truncation 

343  pseudofrequency  s  >  1.  Let  the  first  tail  function  Vq(x)  be 

344  defined  via  (2.27)-(2.29).  Let  cr  E  (0, 1)  be  the  level  of  the  error 

345  in  the  boundary  data,  i.e., 


where  the  number  uj  E  (0, 1)  is  independent  of  n,  k ,  h,  £r  ’"A  355 
and  e*.  356 

Therefore,  Theorem  1  guarantees  that,  if  the  total  number  357 

(n  k) 

Q  of  computed  functions  ;  is  fixed  and  error  parameters  358 
cr,  h  are  sufficiently  small,  then  obtained  iterative  solutions  359 
£rn’^  ( x )  are  sufficiently  close  to  the  exact  solution  e* ,;  and  this  360 
closeness  is  defined  by  the  error  parameters.  Therefore,  the  total  361 
number  of  iterations  Q  can  be  considered  as  the  regularization  362 
parameter  of  our  process,  which  is  the  additional  regularization  363 
parameter  to  the  number  s.  The  combination  of  inequalities  364 
in  (2.30)  and  (2.31)  has  a  direct  analog  in  the  inequality  in  365 
[11,  Lemma  6.2,  p.  156]  for  classical  Landweber  iterations,  366 
which  are  defined  for  a  substantially  different  ill-posed  prob-  367 
lem.  As  to  the  total  number  of  iterations  Q  being  a  regulariza-  368 
tion  parameter  here,  there  is  no  surprise  in  this.  Indeed,  it  is  369 
stated  on  [1 1,  p.  157]  that  the  number  of  iterations  can  serve  as  370 
a  regularization  parameter  for  an  ill-posed  problem.  371 


III.  Imaging  Results 


372 


The  schematic  of  the  data  collection  by  the  forward-looking  373 
radar  is  shown  in  Fig.  2(a).  Time-resolved  electromagnetic  374 
pulses  are  emitted  by  two  sources  installed  on  the  radar.  Only  375 
one  component  of  the  electric  field  is  both  transmitted  and  376 
measured  in  the  backscatter  direction.  The  data  are  collected  377 
by  sixteen  detectors  with  the  step  size  in  time  of  0.133  ns.  378 
Data  from  shallow  targets  placed  both  below  and  above  the  379 
ground  were  provided.  The  only  piece  of  information  provided  380 
by  the  ARL  team  (Sullivan  and  Nguyen)  to  Kuzhuget,  Beilina,  381 
Klibanov,  and  Fiddy  was  whether  the  target  was  located  above  382 
the  ground  or  was  buried.  The  depth  of  the  upper  surface  of  a  383 
buried  target  was  a  few  centimeters.  GPS  was  used  to  provide  384 
the  distance  between  the  radar  and  a  point  on  the  ground,  which  385 
is  located  above  that  target  to  within  a  few  centimeters  error.  386 
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Fig.  3.  (a)  Scattered  field  from  a  metallic  target,  (b)  Scattered  field  from  a  high-permittivity  target  with  the  same  shape  (er  (target)  =  10).  Note  the  similarity 

between  the  backscatter  electric  fields  in  cases  (a)  and  (b). 


387  The  time-resolved  voltages  induced  by  the  backreflected  signals 

388  were  integrated  over  the  radar  to  target  distances  ranging  from  8 

389  to  20  m,  and  they  were  also  averaged  with  respect  to  both  source 

390  positions  and  with  respect  to  the  output  of  the  16  detectors. 

391  Since  we  can  assume  here  that  the  radar/target  distance  was 

392  known,  then  it  was  also  approximately  known  which  part  of  the 

393  measured  time-resolved  signal  would  correspond  to  scattering 

394  events  from  that  target  (see  Fig.  2).  Despite  the  presence  of 

395  clutter,  a  single  time-dependent  curve  is  extracted  from  the 

396  measured  return  time  histories,  as  illustrated  in  Fig.  2(b).  This 

397  is  the  form  of  the  data  that  have  been  processed  in  each  of 

398  the  five  measured  data  sets  generated  by  the  ARL.  A  typical 

399  plot  of  returns  without  applying  the  inverse  algorithm  is  shown 

400  in  Fig.  2(c),  where  the  triangle  denotes  a  possible  target  of 

401  interest  among  the  clutter  from  the  backscatter  generated  from 

402  the  volume  of  the  region  illuminated  by  the  radar  in  Fig.  2(a). 

403  We  process  a  set  of  averaged  time  histories  like  those  shown  in 

404  Fig.  2(b)  to  create  a  down-range  cut  of  the  permittivity  profile, 

405  as  indicated  in  Fig.  2(d). 

406  Our  objective  was  to  calculate  ratios 


R  = 


er  (target) 
n  (background) 


(3.1) 


illustrates  this.  We  use  the  upper  bound  sr  (target)  =  30  for  425 
the  metallic  targets  because  our  calculations  show  that  LT  in  426 
(2.7),  from  the  response  function  g(t),  almost  coincides  for  427 
er  (target)  >  30.  428 

In  both  cases  of  a  metal  structure  and  a  high-permittivity  429 
structure,  one  can  expect  enhanced  backscatter  if  the  incident  430 
pulse  includes  frequencies  that  correspond  to  a  normal  mode  of  431 
the  target.  Hence,  we  assign  432 


10  <  er  (metallic  target)  <  30. 


(3.2) 


407  of  dielectric  constants.  If  the  er  (background)  is  known,  then  it 

408  is  trivial  to  deduce  er  (target) .  Clearly,  for  a  target  located  above 

409  the  ground,  er  (background)  =  1.  In  general,  we  would  expect 

410  the  target  volume  to  contain  many  inhomogeneities  with  spa- 

411  tially  varying  sr(x).  A  weighted  average  of  dielectric  constants 

412  of  these  constituent  materials  will  be  found  over  the  volume 

413  spatial  resolution  cell  that  corresponds  to  the  particular  data 

414  acquisition  configuration.  In  the  examples  presented  here,  we 

415  show  results  obtained  from  just  one  time-history  curve  for  each 

416  target,  corresponding  to  one  polarization  component  of  the  in- 

417  cident  electromagnetic  field  and  backscatter  data  measured  and 

418  averaged  over  all  16  receiver  locations.  Clearly,  this  severely 

419  limits  the  transverse  resolution  but  improves  the  signal-to-noise 

420  ratio  for  1-D  imaging  in  the  propagation  direction.  The  model 

421  is  further  simplified  by  using  the  1-D  CIP  employing  only 

422  one  hyperbolic  PDE.  Consequently,  the  interpretation  of  the 

423  backscattering  radiation  will  assign  a  high-permittivity  value 

424  to  metal  structures.  A  comparison  between  Fig.  3(a)  and  (b) 


We  call  (3.2)  the  appearing  dielectric  constant  of  metallic  tar¬ 
gets.  In  other  words,  we  consider  in  (3.2)  that  regions  appearing 
to  have  a  high  dielectric  constant  could  also  be  metallic  targets. 

To  appreciate  the  kind  of  backscatter  data  and  image  recov¬ 
ery  expected  from  a  simple  dielectric  block,  a  1-D  example 
illustrated  in  Fig.  3  was  investigated.  Computations  in  this 
example  were  performed  using  the  software  package  WavES 
[24].  The  permittivity  profile,  i.e.,  sr  (target)  =  4,  is  shown  in 
Fig.  4(a);  and  the  computed  function  u(0,t)  =  g(t)  for  0  < 
t  <  3  is  shown  in  Fig.  4(b)  [see  (2.3)  for  g(t)].  We  assume 
temporal  units  here  for  which  at  t  =  3,  a  distance  of  x  =  3 
units  is  traversed;  the  source  is  at  xq  =  —  1,  and  the  block’s 
front  face  is  at  x  =  0.2.  Since  the  block  is  0.2  units  wide,  g(t) 
represents  the  backscatter  return  from  the  front  and  back  face  of 
the  block.  The  reason  why,  in  Fig.  4(b),  g(t )  =  0  for  t  <  1  and 
g(t)  =  1/2  for  1  <t<  1.4  is  that  the  solution  of  the  problem  in 
(2.1)  and  (2.2)  for  er(x)  =  1  is  uq(x ,  t)  =  H(t  —  \x  —  x0|)/2, 
where  H(z )  is  the  Heaviside  function,  i.e., 
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440 

441 
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443 

444 

445 

446 

447 

448 

449 

450 


H(z)  = 


(  0,  z  <  0 
\  1,Z>  0. 


Hence,  u(0,t)  =  g(t)  =  H(t  —  l)/2  for  1  <  t  <  1.4;  and  at  451 
t  =  1.4,  the  return  wave  from  the  block  hits  the  observation  452 
point  {x  =  0}  for  the  first  time.  453 

The  measured  data  are  also  challenging  to  process  since  454 
they  arise  from  oblique  illumination,  and  the  exact  location  455 
and  the  amplitudes  of  the  incident  pulses  were  not  known.  In  456 
addition,  a  comparison  of  Fig.  4(b)  with  Fig.  5(b),  (d),  and  (f)  457 
shows  that  the  measured  data  are  highly  oscillatory,  which  are  458 
unlike  their  simulated  counterparts.  Consequently,  we  applied  459 


KUZHUGET  et  al.:  IMAGE  RECOVERY  FROM  BLIND  BACKSCATTERED  DATA  USING  AN  INVERSE  METHOD 


7 


Fig.  4.  (a)  Function  er  (target)  =  4;  note  that  er  (background)  =  1.  (b)  u(0,  t)  =  g(t)  for  0  <  t  <  3.0.  The  source  is  located  at  xq  =  — 1,  and  the  first 

backscatter  return  is  therefore  shown  at  approximately  t  =  2.4  with  “ringing”  determined  by  interference  of  multiply  scattered  waves  between  the  two  boundaries 
of  the  block.  Computations  were  performed  using  the  software  package  WavES  [24]. 
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Fig.  5.  Three  targets  and  their  associated  measured  data.  The  ground  is  dry  sand  with  3  <  er  <  5  [21],  [22].  The  information  shown  in  (a),  (c),  and  (e)  were 
only  provided  after  computations  were  made,  (a)  Depicts  a  bush  that  was  located  on  a  road,  which  generated  background  clutter,  (b)  Scaled  experimental  data  for 
(a),  where  the  horizontal  axis  represents  time  in  nanoseconds  having  a  time  step  of  0.133  ns;  and  the  vertical  axis  is  the  amplitude  of  the  measured  voltage  at  the 
detector,  (c)  Wooden  stake,  (d)  Scaled  experimental  data  for  (c).  (e)  Metal  box  buried  in  dry  sand,  (f)  Scaled  experimental  data  for  (e).  The  mismatch  between 
experimental  and  simulated  data  [see  Fig.  4(b)]  is  evident. 
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Fig.  6.  Superimposed  preprocessed  data  for  all  five  cases  under  consideration. 
The  upward-looking  peak  corresponds  to  the  plastic  cylinder  (see  Table  I). 

460  an  intuitively  reasonable  data  preprocessing  procedure,  which 

461  remained  totally  unbiased  since  it  was  applied  to  blind  data  sets. 

462  The  idea  of  this  procedure  is  to  make  the  data  more  similar  to 

463  that  shown  in  Fig.  4(b).  Previously,  a  similar  procedure  was 

464  reported  for  transmitted  data  in  [6],  [7],  and  [12].  We  have 

465  considered  two  cases. 

466  Case  1.  Suppose  that  the  target  is  located  above  the  ground. 

467  In  this  case 

er (target)  >  sr (background)  =  £r(air)  =  1.  (3.3) 

468  Fig.  4(a)  and  (b)  shows  that,  in  this  case,  the  backscattering 

469  signal  should  be  basically  one  downward-looking  peak. 

470  Therefore,  we  have  selected  on  the  experimental  curve  the 

47 1  first  downward-looking  peak  with  the  largest  amplitude.  As 

472  to  the  rest  of  the  experimental  curve,  it  was  set  to  zero. 

473  Hence,  we  work  only  with  the  selected  peak. 

474  Case  2.  Suppose  that  the  target  is  buried  in  the  ground.  In  this 

475  case,  we  cannot  claim  the  validity  of  (3.3).  On  the  other 

476  hand,  our  numerical  simulations  (not  shown  here)  have 

477  demonstrated  that,  if  er (target)  <  er (background),  then 

478  in  the  analog  of  Fig.  4(b),  the  peak  would  look  upward. 

479  Therefore,  in  this  case,  we  have  selected  on  the  experimen- 

480  tal  curve  of  the  first  peak  with  the  largest  amplitude  to  work 

481  with  initially. 

482  We  were  provided  with  five  data  sets.  Fig.  6  shows  superim- 

483  posed  preprocessed  curves  for  all  five  targets  we  have  worked 

484  with.  The  only  peak  that  looks  upward  is  the  one  for  the  plastic 

485  cylinder  buried  in  soil  since  its  dielectric  constant  was  less 

486  than  that  of  the  soil  (see  Fig.  6).  We  stress  once  again  that 

487  nothing  was  known  in  advance  about  the  dielectric  constants 

488  of  targets.  Therefore,  the  choice  of  the  upward-looking  peak 

489  in  the  case  of  the  plastic  cylinder  was  unbiased  and  was  done 

490  only  using  the  aforementioned  rule.  The  measured  amplitude 

491  for  each  case  was  on  the  order  of  105.  This  is  well  above  the 

492  amplitude  in  Fig.  4(b).  Thus,  all  signals  were  preprocessed 

493  first  (as  previously  described)  and  multiplied  by  the  scaling 


Fig.  7.  Graphs  of  the  function  m(0,  s)  =  w( 0,  s)  —  mo(0,  s)  for  s  E  [1,  5] 
for  the  LT  of  the  computationally  simulated  data  in  Fig.  4(b)  and  of  the 
preprocessed  signal  for  the  bush  (see  Fig.  6).  The  signal  in  Fig.  6  for  the  bush 
was  multiplied  by  10-7.  Minimal  and  maximal  values  of  the  function  w(0:s) 
are  similar  for  both  curves.  A  similar  observation  was  made  for  four  other 
targets  we  have  worked  with. 

number  SN  =  1CT7  next.  Consider  the  LT  of  the  simulated  494 
data  shown  in  Fig.  4(b)  and  then  the  pre-processed  signal  for  the  495 
bush  (see  Fig.  6)  and  multiply  it  by  1CT7.  Fig.  7  depicts  super-  496 
imposed  graphs  of  the  function  w( 0,  s)  =  w(0,  s)  —  ico(0,  s )  497 
for  s  G  [1,  5]  for  both  cases.  One  can  see  that  maximum  and  498 
minimum  values  of  both  curves  are  approximately  the  same.  499 
We  initially  used  SN  =  1(T6,  SN  =  1(T7,  and  SN  =  1(T8.  500 
Only  for  SN  =  10-7,  the  maximum  and  minimum  values  of  501 
functions  w( 0,  s )  for  sG  [1,5]  of  both  aforementioned  curves,  502 
i.e.,  the  one  for  the  LT  of  the  function  depicted  in  Fig.  6  (bush),  503 
being  multiplied  by  10-7,  and  the  one  for  the  LT  of  the  function  504 
in  Fig.  4(b),  were  approximately  the  same.  On  the  other  hand,  505 
those  minimal  and  maximal  values  were  quite  different  from  the  506 
values  of  the  LT  of  the  function  in  Fig.  4(b)  for  SN  =  10-6  and  507 
SN  =  1CT8.  Using  SN  =  10-7,  which  is  based  on  the  data  for  508 
the  bush,  we  have  multiplied  the  other  four  preprocessed  signals  509 
(see  Fig.  6)  by  10-7  and  observed  a  similar  behavior  for  the  four  510 
other  targets.  For  the  case  of  the  inverted  peak  in  Fig.  6,  we  511 
compared  \w(0,  s)|  for  it  with  w( 0,  s )  for  the  aforementioned  512 
simulated  data.  Note  that  the  signals  shown  in  Fig.  6  are  not  yet  513 
multiplied  by  the  scaling  number.  After  multiplying  these  data  514 
by  the  scaling  factor  10-7,  then  for  each  set  of  experimental  515 
data,  we  took  the  resulting  curve  as  the  function  u(0,t)  —  516 
^o(0,t)  :=  g(t)  —uo(0,t).  Next,  we  worked  only  with  this  517 
function  as  the  data,  using  the  aforementioned  algorithm.  For  518 
simple  isolated  targets,  these  steps  of  data  preprocessing  are  519 
justified,  given  the  accuracy  of  the  results  obtained  upon  a  520 
posteriori  inspection.  For  more  complex  target  volumes,  a  more  521 
sophisticated  analysis  of  sets  of  time  histories  will  be  necessary.  522 
The  data  sets  were  processed,  and  the  targets  are  illustrated  in  523 
Fig.  5.  If  we  compare  the  highly  oscillatory  curves  of  Fig.  5(b),  524 
(d)  and  (f),  one  can  see  that  these  backscatter  time  histories  or  525 
signatures  are  qualitatively  quite  similar  in  appearance.  Their  526 
oscillatory  nature  is  due  to  the  specific  carrier  frequency  and  527 
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TABLE  I 

Computed  Values  for  R,  the  Relative  Dielectric  Constant  in  (3.1),  Based  on  Blind  Processing  of  Measured 
Backscatter  Data  From  Five  Different  Targets.  Here,  A  Means  Air,  and  B  Means  Dry  Sand 


Target 

A/B 

R 

er  (backgr) 

er  (target) ,  calc. 

£r  (target) ,  published. 

Figure  3.3-(a) 

n/a 

3.8 

i 

3.8 

4  (known) 

Bush 

A 

6.5 

i 

6.5 

3  to  20  [10] 

Wood  stake 

A 

3.8 

i 

3.8 

2  to  6  [21] 

Metal  box 

B 

3.8 

3  to  5  [21] 

11.4  to  19 

10  to  30  (3.2) 

Metal  cylinder 

B 

4.3 

3  to  5  [21] 

12.9  to  21.4 

10  to  30  (3.2) 

Plastic  cylinder 

B 

0.4 

3  to  5  [21] 

1.2  to  2 

1.1  to  3.2  [21,  22] 

528  finite  bandwidth  of  the  pulsed  radiation,  whereas  the  simulated 

529  data  assume  an  idealized  pulse.  For  these  simple  targets,  we 

530  allow  the  aforementioned  preprocessing  step  to  force  a  cor- 

531  respondence  between  the  two  in  order  to  identify  the  earliest 

532  return  from  the  boundary  of  the  target  and  determine  its  relative 

533  amplitude.  Based  on  this,  the  inversion  algorithm  can  determine 

534  a  reliable  estimate  of  that  target’s  actual  permittivity.  In  addi- 

535  tion,  we  have  conducted  a  limited  sensitivity  study  with  respect 

536  to  the  scaling  factor.  Specifically,  we  took  SN  =  0.8  •  1CT7 

537  and  SN  =  1.2  •  10“ 7  for  all  five  targets,  which  are  variations 

538  of  20%  of  the  scaling  number.  In  five  out  of  five  cases  of 
AQ6  539  experimental  data,  we  have  worked  with  values  of  R  kept  within 

540  tabulated  limits  (see  Table  I)  when  these  variations  of  SN 

541  were  tried.  An  optimal  value  of  SN  might  be  determined  via 

542  a  comparison  of  values  of  R  :=  R(SN)  with  measured  values 

543  for  a  few  known  targets.  At  present,  we  have  concentrated  on 

544  reconstructing  a  real  parameter  that  describes  the  permittivity 

545  of  target  features;  and  metal  objects  have  been  images  simply 

546  having  a  very  large  relative  permittivity.  We  note  that  there  is 

547  no  reason  why  a  conductivity  term  could  not  be  incorporated 

548  into  the  algorithm. 

549  In  addition  to  high  oscillations  of  the  data,  we  have  faced 

550  two  more  uncertainties.  First,  we  did  not  know  where  the 

551  time  t  =  0  is  on  our  data.  Second,  we  did  not  know  where 

552  the  actual  location  of  the  source  xo  is.  This  means  that  it  is 

553  impossible  to  determine  the  location  of  the  target.  Hence,  for 

554  computational  purposes,  we  have  arbitrarily  assigned  t  =  0  to 
AQ7  555  be  a  fixed  location  1  ns  off  to  the  left  from  the  beginning  of  the 

556  largest  amplitude  peak  and  xq  :=  —1,  knowing  that  we  have 

557  independent  GPS  data  to  better  fix  absolute  ranges  should  we 

558  need  that  information.  Our  primary  objective  here  is  to  confirm 

559  the  quantitative  accuracy  of  the  estimates  of  the  dielectric 

560  constant  of  each  of  the  targets,  i.e.,  to  accurately  image  the  ratio 

561  R  in  (3.1). 

562  The  derivative  of  the  LT  of  the  preprocessed  data  was  corn- 

563  puted  for  0  <  s  <  12  with  a  step  size  of  As  =  0.05.  Since 

564  the  calculation  of  the  derivative  of  noisy  data  is  an  ill-posed 

565  problem,  we  have  used  the  following  well-known  formula  for 

566  the  calculation  of  the  derivative  of  the  LT: 

oo 

<p'(s)  -  dsw0(0,  s)  =  -  J  ( g(t )  -  uq (0,  t))  te~stdt.  (3.4) 
o 

567  Since  for  all  targets  the  function  g(t)  —  ^o(0,  t)  =  0  for  t  >  2 

568  (see  Fig.  6),  then  the  integration  in  (3.4)  is  actually  carried  for 

569  0  <  t  <  2.  We  then  define  boundary  conditions  for  functions 

570  qn,k  f°r  each  n,  and  R  is  calculated  by  the  aforementioned 

571  algorithm. 


In  Fig.  8(a)  and  (f),  we  regard  R  as  the  maximal  amplitude  of  572 
the  calculated  peak.  We  first  verified  that  the  algorithm  provides  573 
a  good  estimate  for  R  using  simulated  data.  For  the  block  in  574 
Fig.  4(a),  we  obtain  the  1-D  image  shown  in  Fig.  8(a),  which  575  aqs 
was  found  to  be  er  =  3.8,  which  is  very  close  to  the  known  576 
value  of  4.  Next,  we  have  calculated  images  from  experimental  577 
data.  In  addition  to  Fig.  5(a),  (c),  and  (e),  we  have  also  imaged  578 
two  more  cases,  namely,  a  plastic  cylinder  and  a  metal  cylinder,  579 
which  are  both  buried  in  the  ground  with  schematics  similar  580 
with  the  one  in  Fig.  5(e).  Fig.  8(b)— (f)  displays  our  calculated  581 
images  for  all  five  targets.  582 

Dielectric  constants  were  not  measured  when  the  data  were  583 
collected.  Therefore,  we  have  compared  computed  values  of  584 
dielectric  constants  with  those  listed  in  tables  [21],  [22].  Note  585 
that  these  tables  often  provide  a  range  of  values  rather  than  586 
exact  numbers;  but  given  this  caveat,  the  calculated  results  587 
for  these  materials  are  well  within  the  range  of  expectations  588 
(see  Table  I).  589 

IV.  Conclusion  590 

We  have  described  a  new  method  for  recovering  quanti-  591 
tatively  reliable  estimates  of  target’s  material  properties  (di-  592 
electric  constants)  from  backscatter  field  measurements.  The  593 
method  is  an  inverse  scattering  algorithm  based  on  a  rigorously  594 
formulated  CIP.  The  numerical  method  is  constructed  to  ensure  595 
global  convergence,  and  therefore,  it  avoids  stagnation  at  erro-  596 
neous  solutions  for  images  of  target  permittivity  distributions.  597 
Furthermore,  the  method  requires  no  prior  knowledge  of  the  598 
inhomogeneities  present  in  the  target  volume.  These  properties  599 
are  rigorously  guaranteed.  The  authors  are  unaware  of  alterna-  600 
tive  numerical  methods  with  similar  characteristics  for  the  case  601 
of  the  CIPs  making  use  of  such  limited  data.  602 

The  approach  was  evaluated  here  using  data  provided  by  the  603 
ARL  from  a  forward-looking  radar  system  without  any  prior  604 
knowledge  of  the  targets  being  used.  The  data  were  measured  605 
using  oblique  incidence  and  with  unknown  source  locations,  606 
and  thus,  some  assumptions  were  made  to  provide  the  necessary  607 
inputs  for  the  algorithm.  The  procedure  first  estimates  a  solution  608 
that  has  defined  error  given  the  quality  of  the  data  but  which  609 
is  guaranteed  to  be  reliable.  To  simplify  matters,  only  images  610 
of  dielectric  constants  were  recovered  in  order  to  validate  the  611 
quantitative  accuracy  of  the  approach.  Data  sets  were  prepro-  612 
cessed,  and  a  downrange  permittivity  profile  was  calculated.  613 
If  the  angular  spread  of  backscatter  time  histories  would  be  614 
measured,  then  its  additional  processing  would  provide  a  3-D  615 
image  with  a  high  spatial  resolution,  despite  the  use  here  of  a  616 
single  source  point  (see  [7,  Fig.  6.3]).  617 
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(c)  (d) 


(e)  (f) 


Fig.  8.  Calculated  images  of  targets.  The  ratio  R  in  (3.1)  is  regarded  as  the  maximal  amplitude  of  the  imaged  peak,  (a)  Image  for  computationally  simulated  data 
as  a  verification  of  the  accuracy  of  our  algorithm.  The  rectangular  block  and  the  curve  are  true  and  computed  profiles  of  the  dielectric  constant,  respectively.  The 
computed  target/background  contrast  R  =  3.8,  which  corresponds  to  a  5%  of  error,  (b)  Image  of  the  bush  [see  Fig.  2(a)].  The  calculated  er  (bush)  =  6.5,  which 
is  in  the  range  of  tabulated  values  3  <  er  <  20  [10].  (c)  Image  of  the  wood  stake  [see  Fig.  4(c)].  The  calculated  er(wood  stake)  =  3.8  [10].  (d)  Image  of  the 
buried  metal  box  [see  Fig.  5(e)].  The  calculated  R  =  3.8.  Since  the  background  was  dry  sand  with  3  <  er  (dry  sand)  <  5  [21],  then  the  calculated  er  (metal  box) 
is  between  11.4  and  19.  This  is  within  the  range  [see  (3.2)]  of  appearing  dielectric  constants  of  metallic  targets,  (e)  Calculated  image  of  the  buried  metal  cylinder. 
The  calculated  ratio  R  =  4.3.  Similarly  with  (d),  we  conclude  that  the  calculated  value  of  er  (metal  cylinder)  is  between  12.9  and  21.4.  This  is  again  within  the 
range  [see  (3.2)]  of  appearing  dielectric  constants  of  metallic  targets,  (f)  Calculated  image  of  the  buried  plastic  cylinder.  The  calculated  ratio  R  =  0.4.  Similarly 
with  (d),  we  conclude  that  the  calculated  value  of  the  dielectric  constant  er  (plastic  cylinder)  is  between  1.2  and  2.5,  which  is  again  within  the  range  of  tabulated 
values  for  plastic  [21],  [22]. 
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618  Since  the  dielectric  constants  of  the  targets  were  not  actually 

619  measured  in  the  ARL  experiments,  then  the  best  one  can  do 

620  is  compare  retrieved  parameters  with  tabulated  values.  Table  I 

621  shows  the  computed  relative  permittivities  of  targets.  It  is 

622  clearly  shown  that  all  five  targets  fall  well  within  expected 

623  tabulated  limits  for  the  materials  in  question,  despite  the  fact 

624  that  no  prior  knowledge  whatsoever  was  employed.  We  further 

625  emphasize  that  these  results  were  obtained  despite  a  very  lim- 

626  ited  information  content,  large  noise  in  the  data,  and  significant 

627  discrepancies  between  experimental  and  simulated  data.  We  can 

628  therefore  conclude  that  these  results  point  toward  the  validity  of 

629  our  mathematical  model.  The  fact  that  regardless  of  limitations 

630  of  the  method,  we  consistently  got  results,  which  only  later 

631  were  found  to  fall  well  within  tabulated  limits,  points  toward 

632  a  great  degree  of  robustness  of  this  algorithm. 

633  The  purpose  of  estimating  the  dielectric  constant  is  to  provide 

634  one  extra  piece  of  information  about  the  target.  Up  to  this  point, 

635  most  of  the  radar  community  has  solely  relied  on  the  intensity 

636  of  the  radar  image  for  doing  detection  and  discrimination.  It  is 

637  anticipated  that,  when  the  intensity  information  is  coupled  with 

638  the  new  dielectric  information  this  method  provides,  algorithms 

639  can  be  then  designed  that  will  provide  better  performance  in 

640  terms  of  probability  of  detection  and  false  alarm  rates.  Finally, 

641  we  repeat  that  the  results  presented  in  this  paper  are  primarily 

642  being  used  as  a  vehicle  to  illustrate  this  powerful  inverse 

643  scattering  algorithm  method  and  its  ability  to  recover  dielectric 

644  properties  of  targets  from  experimental  data  collected  by  the 

645  forward-looking  radar  of  the  ARL.  Detailed  studies  making 

646  use  of  larger  experimental  data  sets  from  more  complex  3-D 

647  scattering  objects  are  necessary,  and  the  authors  will  report  on 

648  this  in  the  near  future. 
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